CovidFinal
Daniel Clark
7/29/2020
File and Library Load
library(tswge)
# Read in the Arizona Data
azData = read.csv('/Users/danielclark/Desktop/SMU/Time_Series/ImportsProject/Project/Github/6373TimeSeries_Covid-master/original_data/AZdaily.csv', header=TRUE)
# Read in the United States Data
usData = read.csv('/Users/danielclark/Desktop/SMU/Time_Series/ImportsProject/Project/Github/6373TimeSeries_Covid-master/original_data/USdaily.csv', header=TRUE)Data Cleaning
# Sort the data by date
as2 = azData[order(azData$date),]
us2 = usData[order(usData$date),]
head(as2)## date state positive negative pending hospitalizedCurrently
## 140 20200304 AZ 2 27 5 NA
## 139 20200305 AZ 2 28 6 NA
## 138 20200306 AZ 3 33 15 NA
## 137 20200307 AZ 5 44 7 NA
## 136 20200308 AZ 5 44 7 NA
## 135 20200309 AZ 5 44 7 NA
## hospitalizedCumulative inIcuCurrently inIcuCumulative onVentilatorCurrently
## 140 4 NA NA NA
## 139 5 NA NA NA
## 138 6 NA NA NA
## 137 6 NA NA NA
## 136 6 NA NA NA
## 135 8 NA NA NA
## onVentilatorCumulative recovered dataQualityGrade lastUpdateEt dateModified
## 140 NA NA
## 139 NA NA
## 138 NA NA
## 137 NA NA
## 136 NA NA
## 135 NA NA
## checkTimeEt death hospitalized dateChecked totalTestsViral
## 140 NA 4 29
## 139 NA 5 30
## 138 NA 6 36
## 137 NA 6 49
## 136 NA 6 49
## 135 NA 8 49
## positiveTestsViral negativeTestsViral positiveCasesViral deathConfirmed
## 140 NA NA NA NA
## 139 NA NA NA NA
## 138 NA NA NA NA
## 137 NA NA NA NA
## 136 NA NA NA NA
## 135 NA NA NA NA
## deathProbable fips positiveIncrease negativeIncrease total totalTestResults
## 140 NA 4 0 0 34 29
## 139 NA 4 0 1 36 30
## 138 NA 4 1 5 51 36
## 137 NA 4 2 11 56 49
## 136 NA 4 0 0 56 49
## 135 NA 4 0 0 56 49
## totalTestResultsIncrease posNeg deathIncrease hospitalizedIncrease
## 140 0 29 0 0
## 139 1 30 0 1
## 138 6 36 0 1
## 137 13 49 0 0
## 136 0 49 0 0
## 135 0 49 0 2
## hash commercialScore
## 140 7b439d158b4705a4482d78078c5fbf7349d455b9 0
## 139 c34463e6bdd1606ff5a91b3a2fb86f7d115ee32c 0
## 138 cedf256abd5258aeb7d8f878405935ebc92d3501 0
## 137 a9cdd6464832a95ece09ed796584dc5c80e859d0 0
## 136 fc5d4224c90b1a71a22d0c5cad290490cc2528be 0
## 135 0061536552e197c87f2b29a428cd4d44ea0fb1b2 0
## negativeRegularScore negativeScore positiveScore score grade
## 140 0 0 0 0 NA
## 139 0 0 0 0 NA
## 138 0 0 0 0 NA
## 137 0 0 0 0 NA
## 136 0 0 0 0 NA
## 135 0 0 0 0 NAhead(us2)## date states positive negative pending hospitalizedCurrently
## 182 20200122 1 2 0 NA NA
## 181 20200123 1 2 0 NA NA
## 180 20200124 1 2 0 NA NA
## 179 20200125 1 2 0 NA NA
## 178 20200126 1 2 0 NA NA
## 177 20200127 1 2 0 NA NA
## hospitalizedCumulative inIcuCurrently inIcuCumulative onVentilatorCurrently
## 182 NA NA NA NA
## 181 NA NA NA NA
## 180 NA NA NA NA
## 179 NA NA NA NA
## 178 NA NA NA NA
## 177 NA NA NA NA
## onVentilatorCumulative recovered dateChecked death hospitalized
## 182 NA NA 2020-01-22T00:00:00Z NA NA
## 181 NA NA 2020-01-23T00:00:00Z NA NA
## 180 NA NA 2020-01-24T00:00:00Z NA NA
## 179 NA NA 2020-01-25T00:00:00Z NA NA
## 178 NA NA 2020-01-26T00:00:00Z NA NA
## 177 NA NA 2020-01-27T00:00:00Z NA NA
## lastModified total totalTestResults posNeg deathIncrease
## 182 2020-01-22T00:00:00Z 2 2 2 0
## 181 2020-01-23T00:00:00Z 2 2 2 0
## 180 2020-01-24T00:00:00Z 2 2 2 0
## 179 2020-01-25T00:00:00Z 2 2 2 0
## 178 2020-01-26T00:00:00Z 2 2 2 0
## 177 2020-01-27T00:00:00Z 2 2 2 0
## hospitalizedIncrease negativeIncrease positiveIncrease
## 182 0 0 0
## 181 0 0 0
## 180 0 0 0
## 179 0 0 0
## 178 0 0 0
## 177 0 0 0
## totalTestResultsIncrease hash
## 182 0 d538c99729d1fee626212d1878a100c1e1204a5f
## 181 0 cee36ebf3174bf1df0daa36e1e8088a157406fad
## 180 0 bfffe76fc0b7cf11efe8aecd3cc7b22598d77d61
## 179 0 bef2a1d5f2a13491e0e0369bbd46c10cdd12973b
## 178 0 e1cf59ab48e1cf367c4a6798a508a23d9d36bd18
## 177 0 3dde2425e73c17b4e6cbf0e4432340298c0a94f5Plot the US data to make sure it all looks good.
# Plot of date, just as a QA/QC Check - all looks good.
plotts.sample.wge(us2$date)
## $autplt
## [1] 1.0000000 0.9812540 0.9624335 0.9435387 0.9245698 0.9055269 0.8864103
## [8] 0.8672200 0.8479564 0.8286195 0.8092096 0.7935956 0.7779212 0.7621866
## [15] 0.7463920 0.7305375 0.7146234 0.6986498 0.6826170 0.6665249 0.6503740
## [22] 0.6341643 0.6214856 0.6087610 0.5959906 0.5831747
##
## $freq
## [1] 0.005494505 0.010989011 0.016483516 0.021978022 0.027472527 0.032967033
## [7] 0.038461538 0.043956044 0.049450549 0.054945055 0.060439560 0.065934066
## [13] 0.071428571 0.076923077 0.082417582 0.087912088 0.093406593 0.098901099
## [19] 0.104395604 0.109890110 0.115384615 0.120879121 0.126373626 0.131868132
## [25] 0.137362637 0.142857143 0.148351648 0.153846154 0.159340659 0.164835165
## [31] 0.170329670 0.175824176 0.181318681 0.186813187 0.192307692 0.197802198
## [37] 0.203296703 0.208791209 0.214285714 0.219780220 0.225274725 0.230769231
## [43] 0.236263736 0.241758242 0.247252747 0.252747253 0.258241758 0.263736264
## [49] 0.269230769 0.274725275 0.280219780 0.285714286 0.291208791 0.296703297
## [55] 0.302197802 0.307692308 0.313186813 0.318681319 0.324175824 0.329670330
## [61] 0.335164835 0.340659341 0.346153846 0.351648352 0.357142857 0.362637363
## [67] 0.368131868 0.373626374 0.379120879 0.384615385 0.390109890 0.395604396
## [73] 0.401098901 0.406593407 0.412087912 0.417582418 0.423076923 0.428571429
## [79] 0.434065934 0.439560440 0.445054945 0.450549451 0.456043956 0.461538462
## [85] 0.467032967 0.472527473 0.478021978 0.483516484 0.489010989 0.494505495
## [91] 0.500000000
##
## $db
## [1] 17.2944492 11.3185302 7.9845353 5.3591595 3.4235257 4.7228467
## [7] 0.6616461 -0.5886179 -1.6146724 -2.7779072 -3.8740248 -0.2254521
## [13] -4.8774155 -5.6853788 -7.1839471 -6.8548119 -6.9696806 -13.4437757
## [19] -8.6136800 -8.6135467 -8.1768315 -8.7801547 -8.9355394 -9.0139596
## [25] -9.9458302 -10.4281262 -9.8246180 -11.3599250 -12.6088499 -7.6473271
## [31] -11.7371524 -12.4102651 -14.8244061 -13.6523576 -13.2471424 -14.8916965
## [37] -14.2182261 -13.7682588 -13.4583520 -13.1824947 -12.5724209 -16.0995791
## [43] -14.3422559 -14.6609758 -12.5841529 -14.4930707 -15.9471203 -11.6332019
## [49] -14.8215164 -15.5470337 -17.5101227 -17.4025248 -17.3669351 -15.4851192
## [55] -16.6989938 -15.9764170 -17.3145281 -15.7662136 -14.3586050 -20.6333307
## [61] -16.8671578 -17.1069073 -14.0619680 -15.6907929 -16.5986587 -14.5989336
## [67] -16.6711901 -17.6307343 -17.4615105 -19.1157865 -20.7142327 -15.9853302
## [73] -17.6686506 -16.5821976 -20.6247623 -17.6541107 -15.5583533 -21.1210527
## [79] -17.9546975 -17.8725588 -15.1523188 -16.0771161 -16.0624623 -16.8392128
## [85] -17.7713005 -19.4906065 -16.5879219 -18.9540186 -22.2946389 -16.4191415
## [91] -17.9094268
##
## $dbz
## [1] 12.2494880 11.8587083 11.2034088 10.2781447 9.0768427 7.5956457
## [7] 5.8392120 3.8338270 1.6512527 -0.5614206 -2.5853404 -4.2343819
## [13] -5.4986069 -6.5113919 -7.3752335 -8.0803123 -8.5710912 -8.8609666
## [19] -9.0609770 -9.3058222 -9.6800060 -10.1946812 -10.7969207 -11.3974837
## [25] -11.9164912 -12.3273719 -12.6610268 -12.9670950 -13.2736234 -13.5764749
## [31] -13.8571823 -14.1097952 -14.3537171 -14.6226182 -14.9404889 -15.3027576
## [37] -15.6730537 -15.9986372 -16.2388836 -16.3879353 -16.4722765 -16.5276116
## [43] -16.5778934 -16.6310867 -16.6889830 -16.7595226 -16.8611250 -17.0163028
## [49] -17.2398092 -17.5286050 -17.8584435 -18.1892811 -18.4789175 -18.6988507
## [55] -18.8420729 -18.9177451 -18.9394827 -18.9183165 -18.8643487 -18.7922729
## [61] -18.7235933 -18.6828035 -18.6903705 -18.7570409 -18.8819871 -19.0547470
## [67] -19.2593495 -19.4783836 -19.6950971 -19.8931901 -20.0561105 -20.1684985
## [73] -20.2205969 -20.2130210 -20.1577401 -20.0737754 -19.9804588 -19.8925845
## [79] -19.8196226 -19.7680577 -19.7441786 -19.7548868 -19.8057552 -19.8974756
## [85] -20.0228907 -20.1666940 -20.3087817 -20.4302739 -20.5190659 -20.5713661
## [91] -20.5884010# Examine some plots of the data, such as the cumulative death rate.
plotts.sample.wge(us2$death[36:182])
## $autplt
## [1] 1.0000000 0.9866492 0.9729699 0.9588458 0.9443278 0.9294891 0.9143373
## [8] 0.8988644 0.8830399 0.8668535 0.8502461 0.8332646 0.8159763 0.7983944
## [15] 0.7805131 0.7623341 0.7438786 0.7250291 0.7057989 0.6862189 0.6663467
## [22] 0.6461943 0.6257629 0.6050456 0.5840129 0.5626699
##
## $freq
## [1] 0.006802721 0.013605442 0.020408163 0.027210884 0.034013605 0.040816327
## [7] 0.047619048 0.054421769 0.061224490 0.068027211 0.074829932 0.081632653
## [13] 0.088435374 0.095238095 0.102040816 0.108843537 0.115646259 0.122448980
## [19] 0.129251701 0.136054422 0.142857143 0.149659864 0.156462585 0.163265306
## [25] 0.170068027 0.176870748 0.183673469 0.190476190 0.197278912 0.204081633
## [31] 0.210884354 0.217687075 0.224489796 0.231292517 0.238095238 0.244897959
## [37] 0.251700680 0.258503401 0.265306122 0.272108844 0.278911565 0.285714286
## [43] 0.292517007 0.299319728 0.306122449 0.312925170 0.319727891 0.326530612
## [49] 0.333333333 0.340136054 0.346938776 0.353741497 0.360544218 0.367346939
## [55] 0.374149660 0.380952381 0.387755102 0.394557823 0.401360544 0.408163265
## [61] 0.414965986 0.421768707 0.428571429 0.435374150 0.442176871 0.448979592
## [67] 0.455782313 0.462585034 0.469387755 0.476190476 0.482993197 0.489795918
## [73] 0.496598639
##
## $db
## [1] 17.4166582 8.7695293 4.2804573 2.3223206 0.9142909 -1.0918255
## [7] -2.2496570 -3.2875079 -4.3627783 -5.3632001 -6.3900921 -6.9492473
## [13] -7.4135876 -8.2179617 -8.6354960 -9.3335360 -9.9432768 -10.0228209
## [19] -10.7852442 -10.5956736 -12.6428284 -11.7905503 -12.5392450 -12.6008901
## [25] -12.8916251 -13.2010364 -13.8361817 -13.9237888 -14.0558054 -14.4070371
## [31] -14.4884995 -14.9950913 -15.1873394 -15.1964089 -15.4206044 -15.5554435
## [37] -15.9562751 -16.1824591 -16.2667044 -16.4977961 -16.3440190 -16.8046610
## [43] -16.8975457 -16.9022516 -17.1429496 -17.0660244 -17.3000060 -17.4305677
## [49] -17.6909851 -17.7086947 -17.7383023 -17.6861594 -17.9656942 -18.1376415
## [55] -18.1448033 -18.1908361 -18.0639892 -18.3886535 -18.2880956 -18.6829472
## [61] -18.4940155 -18.4812072 -18.2852802 -18.7373806 -18.6397754 -18.8528607
## [67] -18.5520486 -18.7205487 -18.6528237 -18.9274114 -18.8751583 -18.8437393
## [73] -18.6581417
##
## $dbz
## [1] 12.0058438 11.4812785 10.5972948 9.3396063 7.6896756 5.6293747
## [7] 3.1554321 0.3195886 -2.6811858 -5.4114980 -7.4081915 -8.6837790
## [13] -9.5527603 -10.1176975 -10.3407515 -10.3446319 -10.3962327 -10.7059484
## [19] -11.3457627 -12.2579451 -13.2774787 -14.1937395 -14.8787034 -15.3599012
## [25] -15.7259759 -16.0062162 -16.1736850 -16.2330006 -16.2655614 -16.3810369
## [31] -16.6457149 -17.0500610 -17.5185753 -17.9540384 -18.2977375 -18.5556656
## [37] -18.7636266 -18.9374024 -19.0636228 -19.1324199 -19.1689126 -19.2262791
## [43] -19.3500012 -19.5486563 -19.7907996 -20.0269738 -20.2220056 -20.3722265
## [49] -20.4934786 -20.5963330 -20.6758593 -20.7239012 -20.7477962 -20.7737753
## [55] -20.8308720 -20.9308913 -21.0608228 -21.1923650 -21.3011440 -21.3804101
## [61] -21.4381738 -21.4830405 -21.5144941 -21.5267390 -21.5206115 -21.5099964
## [67] -21.5151659 -21.5488606 -21.6071025 -21.6718943 -21.7236830 -21.7536210
## [73] -21.7652567# Plot the total number of positive cases, negative cases, and the ratio between he two.
plotts.sample.wge(us2$positive)
## $autplt
## [1] 1.0000000 0.9802808 0.9606923 0.9411760 0.9218403 0.9026998 0.8839126
## [8] 0.8654128 0.8471438 0.8290888 0.8112084 0.7935475 0.7761431 0.7590512
## [15] 0.7421863 0.7256001 0.7091790 0.6928110 0.6765577 0.6605387 0.6447581
## [22] 0.6292205 0.6138976 0.5987622 0.5837244 0.5688268
##
## $freq
## [1] 0.005494505 0.010989011 0.016483516 0.021978022 0.027472527 0.032967033
## [7] 0.038461538 0.043956044 0.049450549 0.054945055 0.060439560 0.065934066
## [13] 0.071428571 0.076923077 0.082417582 0.087912088 0.093406593 0.098901099
## [19] 0.104395604 0.109890110 0.115384615 0.120879121 0.126373626 0.131868132
## [25] 0.137362637 0.142857143 0.148351648 0.153846154 0.159340659 0.164835165
## [31] 0.170329670 0.175824176 0.181318681 0.186813187 0.192307692 0.197802198
## [37] 0.203296703 0.208791209 0.214285714 0.219780220 0.225274725 0.230769231
## [43] 0.236263736 0.241758242 0.247252747 0.252747253 0.258241758 0.263736264
## [49] 0.269230769 0.274725275 0.280219780 0.285714286 0.291208791 0.296703297
## [55] 0.302197802 0.307692308 0.313186813 0.318681319 0.324175824 0.329670330
## [61] 0.335164835 0.340659341 0.346153846 0.351648352 0.357142857 0.362637363
## [67] 0.368131868 0.373626374 0.379120879 0.384615385 0.390109890 0.395604396
## [73] 0.401098901 0.406593407 0.412087912 0.417582418 0.423076923 0.428571429
## [79] 0.434065934 0.439560440 0.445054945 0.450549451 0.456043956 0.461538462
## [85] 0.467032967 0.472527473 0.478021978 0.483516484 0.489010989 0.494505495
## [91] 0.500000000
##
## $db
## [1] 17.7632249 9.5344132 7.8680395 5.7812048 3.4364161 2.0248643
## [7] 0.5823541 -0.6316365 -1.4899993 -2.4258452 -3.2259209 -3.9874132
## [13] -4.6323020 -5.3897323 -5.9365950 -6.4957787 -7.0126163 -7.5231456
## [19] -7.9113257 -8.3359654 -8.7427317 -9.0580077 -9.4548722 -9.8211688
## [25] -10.0766269 -10.5607798 -11.0888054 -11.2762083 -11.5097908 -11.7081175
## [31] -12.0001530 -12.2041428 -12.4844143 -12.6649415 -12.9248617 -13.1208810
## [37] -13.3602059 -13.5314362 -13.7728763 -13.9207899 -14.0853717 -14.2235032
## [43] -14.4343725 -14.5329741 -14.7266274 -14.8322199 -15.0024590 -15.1523735
## [49] -15.2629208 -15.4871294 -15.6273129 -15.6754536 -15.7805922 -15.8331769
## [55] -15.9413121 -16.1213286 -16.1702148 -16.3075178 -16.3774717 -16.5418706
## [61] -16.5723456 -16.7267639 -16.7713714 -16.8341397 -16.9398469 -16.9972700
## [67] -17.0329028 -17.1002365 -17.1574101 -17.2428758 -17.2835378 -17.3462392
## [73] -17.3768339 -17.4557004 -17.5000172 -17.5351932 -17.5719268 -17.5655341
## [79] -17.6225219 -17.6339382 -17.6252487 -17.7075151 -17.6765159 -17.7012743
## [85] -17.7465839 -17.7026476 -17.7694889 -17.7255590 -17.7882101 -17.8170126
## [91] -17.8161472
##
## $dbz
## [1] 12.2411168 11.8507688 11.1961269 10.2716195 9.0709137 7.5896440
## [7] 5.8315204 3.8212438 1.6286174 -0.5993725 -2.6364163 -4.2782973
## [13] -5.4976863 -6.4239824 -7.1681263 -7.7417663 -8.1228103 -8.3503566
## [19] -8.5357524 -8.7986418 -9.2114627 -9.7828804 -10.4638793 -11.1674446
## [25] -11.8053777 -12.3312583 -12.7546548 -13.1114989 -13.4235929 -13.6852560
## [31] -13.8822741 -14.0206098 -14.1363770 -14.2791245 -14.4857372 -14.7638021
## [37] -15.0903399 -15.4241956 -15.7268800 -15.9806734 -16.1913253 -16.3746144
## [43] -16.5401606 -16.6858718 -16.8053155 -16.8998300 -16.9841419 -17.0806291
## [49] -17.2072456 -17.3679373 -17.5510196 -17.7358767 -17.9040563 -18.0479221
## [55] -18.1712266 -18.2821951 -18.3855126 -18.4794868 -18.5597589 -18.6257859
## [61] -18.6843981 -18.7473814 -18.8249743 -18.9200141 -19.0264761 -19.1330833
## [67] -19.2295901 -19.3116470 -19.3811457 -19.4424033 -19.4976525 -19.5454040
## [73] -19.5826604 -19.6089035 -19.6284437 -19.6490823 -19.6779840 -19.7176691
## [79] -19.7647520 -19.8121255 -19.8530986 -19.8847515 -19.9083478 -19.9268655
## [85] -19.9418388 -19.9520291 -19.9548687 -19.9494409 -19.9385386 -19.9280107
## [91] -19.9236925plotts.sample.wge(us2$negative)
## $autplt
## [1] 1.0000000 0.9787712 0.9576538 0.9366386 0.9157710 0.8950489 0.8745720
## [8] 0.8543257 0.8342538 0.8143670 0.7946362 0.7750748 0.7555821 0.7363794
## [15] 0.7173393 0.6983833 0.6795372 0.6606898 0.6420026 0.6234559 0.6051504
## [22] 0.5870282 0.5690659 0.5512769 0.5336780 0.5162231
##
## $freq
## [1] 0.005494505 0.010989011 0.016483516 0.021978022 0.027472527 0.032967033
## [7] 0.038461538 0.043956044 0.049450549 0.054945055 0.060439560 0.065934066
## [13] 0.071428571 0.076923077 0.082417582 0.087912088 0.093406593 0.098901099
## [19] 0.104395604 0.109890110 0.115384615 0.120879121 0.126373626 0.131868132
## [25] 0.137362637 0.142857143 0.148351648 0.153846154 0.159340659 0.164835165
## [31] 0.170329670 0.175824176 0.181318681 0.186813187 0.192307692 0.197802198
## [37] 0.203296703 0.208791209 0.214285714 0.219780220 0.225274725 0.230769231
## [43] 0.236263736 0.241758242 0.247252747 0.252747253 0.258241758 0.263736264
## [49] 0.269230769 0.274725275 0.280219780 0.285714286 0.291208791 0.296703297
## [55] 0.302197802 0.307692308 0.313186813 0.318681319 0.324175824 0.329670330
## [61] 0.335164835 0.340659341 0.346153846 0.351648352 0.357142857 0.362637363
## [67] 0.368131868 0.373626374 0.379120879 0.384615385 0.390109890 0.395604396
## [73] 0.401098901 0.406593407 0.412087912 0.417582418 0.423076923 0.428571429
## [79] 0.434065934 0.439560440 0.445054945 0.450549451 0.456043956 0.461538462
## [85] 0.467032967 0.472527473 0.478021978 0.483516484 0.489010989 0.494505495
## [91] 0.500000000
##
## $db
## [1] 17.4825152 11.2511134 7.6702075 5.5378117 3.3947266 1.8709734
## [7] 0.6188099 -0.7321948 -1.5781578 -2.4786055 -3.2635353 -4.0348195
## [13] -4.5796085 -5.3996115 -6.0040718 -6.4994177 -7.0504426 -7.5894163
## [19] -7.9143949 -8.4429318 -8.7165383 -9.1882333 -9.5855305 -9.8938189
## [25] -10.0435837 -10.5417900 -11.0557728 -11.3334835 -11.4398893 -11.7406284
## [31] -12.0063819 -12.2533423 -12.4869953 -12.6870309 -12.8856104 -13.0689437
## [37] -13.3313633 -13.5428188 -13.6362577 -14.0539614 -14.1357028 -14.2455382
## [43] -14.4148629 -14.6144338 -14.6929502 -14.9274889 -15.0707300 -15.1766953
## [49] -15.2058897 -15.5060006 -15.5233842 -15.5891296 -15.7667279 -15.9075334
## [55] -15.9991479 -16.2154541 -16.1602234 -16.3971541 -16.3912602 -16.6481605
## [61] -16.5608955 -16.7307834 -16.8131594 -16.8876663 -16.8613470 -16.9914401
## [67] -17.0030301 -17.1692692 -17.1085926 -17.2372617 -17.2348685 -17.4261193
## [73] -17.4373428 -17.5469076 -17.5208470 -17.5948440 -17.5848075 -17.6231736
## [79] -17.5512996 -17.6258984 -17.6191894 -17.7564142 -17.6438485 -17.7121090
## [85] -17.7467224 -17.8602215 -17.8117072 -17.8824435 -17.7892197 -17.7933063
## [91] -17.8290714
##
## $dbz
## [1] 12.1935643 11.8090537 11.1645499 10.2551374 9.0754769 7.6225600
## [7] 5.9016893 3.9387309 1.8022160 -0.3684802 -2.3626642 -3.9878814
## [13] -5.2076546 -6.1328994 -6.8709427 -7.4457310 -7.8465823 -8.1086705
## [19] -8.3279170 -8.6106306 -9.0239667 -9.5779687 -10.2296860 -10.9014699
## [25] -11.5143555 -12.0255753 -12.4411414 -12.7913930 -13.0965208 -13.3541673
## [31] -13.5543397 -13.7031439 -13.8316872 -13.9828526 -14.1892871 -14.4585879
## [37] -14.7717138 -15.0934974 -15.3903588 -15.6456087 -15.8621066 -16.0516249
## [43] -16.2215689 -16.3701771 -16.4924948 -16.5904546 -16.6775967 -16.7742040
## [49] -16.8969198 -17.0503437 -17.2254213 -17.4049658 -17.5729091 -17.7215984
## [55] -17.8525748 -17.9712233 -18.0803657 -18.1779782 -18.2603085 -18.3272901
## [61] -18.3854485 -18.4457135 -18.5177443 -18.6048268 -18.7025634 -18.8019043
## [67] -18.8944171 -18.9763323 -19.0488328 -19.1149155 -19.1756808 -19.2290288
## [73] -19.2716287 -19.3024468 -19.3249177 -19.3459287 -19.3722951 -19.4071997
## [79] -19.4488929 -19.4922907 -19.5322364 -19.5661786 -19.5945419 -19.6188496
## [85] -19.6393557 -19.6542024 -19.6608338 -19.6585779 -19.6502904 -19.6415763
## [91] -19.6379171# make a plot of the ratio of negative tests to positive tests.
plotts.sample.wge(us2$negative/us2$positive)
## $autplt
## [1] 1.0000000 0.9880294 0.9754470 0.9619649 0.9475554 0.9327023 0.9170225
## [8] 0.9007174 0.8837430 0.8661743 0.8479350 0.8293817 0.8105104 0.7910318
## [15] 0.7715514 0.7516197 0.7316136 0.7115800 0.6915333 0.6715523 0.6515066
## [22] 0.6313961 0.6112190 0.5909700 0.5707745 0.5505423
##
## $freq
## [1] 0.005494505 0.010989011 0.016483516 0.021978022 0.027472527 0.032967033
## [7] 0.038461538 0.043956044 0.049450549 0.054945055 0.060439560 0.065934066
## [13] 0.071428571 0.076923077 0.082417582 0.087912088 0.093406593 0.098901099
## [19] 0.104395604 0.109890110 0.115384615 0.120879121 0.126373626 0.131868132
## [25] 0.137362637 0.142857143 0.148351648 0.153846154 0.159340659 0.164835165
## [31] 0.170329670 0.175824176 0.181318681 0.186813187 0.192307692 0.197802198
## [37] 0.203296703 0.208791209 0.214285714 0.219780220 0.225274725 0.230769231
## [43] 0.236263736 0.241758242 0.247252747 0.252747253 0.258241758 0.263736264
## [49] 0.269230769 0.274725275 0.280219780 0.285714286 0.291208791 0.296703297
## [55] 0.302197802 0.307692308 0.313186813 0.318681319 0.324175824 0.329670330
## [61] 0.335164835 0.340659341 0.346153846 0.351648352 0.357142857 0.362637363
## [67] 0.368131868 0.373626374 0.379120879 0.384615385 0.390109890 0.395604396
## [73] 0.401098901 0.406593407 0.412087912 0.417582418 0.423076923 0.428571429
## [79] 0.434065934 0.439560440 0.445054945 0.450549451 0.456043956 0.461538462
## [85] 0.467032967 0.472527473 0.478021978 0.483516484 0.489010989 0.494505495
## [91] 0.500000000
##
## $db
## [1] 17.066970 14.037677 6.927412 -2.316635 2.942087 2.537177
## [7] -1.868969 -4.297980 -2.582069 -3.921793 -5.472261 -5.833139
## [13] -6.035251 -8.594486 -7.349269 -6.716720 -9.862472 -10.764354
## [19] -8.666251 -9.493776 -10.565484 -12.899658 -11.172070 -10.276316
## [25] -11.896496 -14.208790 -12.200390 -11.615490 -13.277167 -14.795232
## [31] -13.438488 -13.494326 -13.987071 -15.185772 -14.939137 -14.014920
## [37] -15.325682 -15.584789 -14.663605 -16.943673 -15.686198 -15.037756
## [43] -16.671816 -17.711030 -15.386260 -15.958247 -18.239703 -17.221042
## [49] -15.904360 -17.632213 -17.593731 -17.052220 -17.315252 -18.339674
## [55] -17.818034 -17.228004 -18.366872 -19.514684 -16.606651 -17.811078
## [61] -20.529366 -18.124087 -16.460551 -19.723425 -19.957292 -17.096739
## [67] -18.316973 -20.721487 -18.009053 -18.686544 -19.985315 -18.808006
## [73] -17.862535 -21.233827 -19.681351 -17.041450 -19.319917 -22.191100
## [79] -18.074007 -18.132394 -20.218547 -20.640788 -18.688699 -19.367210
## [85] -19.642057 -19.414524 -20.336997 -20.284099 -18.121614 -19.357137
## [91] -21.352540
##
## $dbz
## [1] 12.360881 11.966699 11.304483 10.366484 9.142323 7.620246
## [7] 5.790508 3.653875 1.241416 -1.345822 -3.885746 -6.058459
## [13] -7.657594 -8.757697 -9.525040 -10.010433 -10.213146 -10.224757
## [19] -10.224360 -10.374610 -10.762430 -11.395586 -12.212589 -13.095667
## [25] -13.906727 -14.553869 -15.036149 -15.409838 -15.717018 -15.955016
## [31] -16.103926 -16.173626 -16.219795 -16.316626 -16.518220 -16.836559
## [37] -17.239441 -17.664771 -18.048948 -18.357486 -18.594863 -18.785240
## [43] -18.944928 -19.071866 -19.157419 -19.206310 -19.245054 -19.311405
## [49] -19.434825 -19.621961 -19.853779 -20.094554 -20.308820 -20.477186
## [55] -20.600266 -20.688978 -20.751469 -20.788230 -20.797955 -20.787153
## [61] -20.773951 -20.782359 -20.831503 -20.927059 -21.059249 -21.207756
## [67] -21.350886 -21.473924 -21.571777 -21.645162 -21.694865 -21.719679
## [73] -21.719727 -21.701682 -21.680512 -21.675354 -21.701984 -21.766417
## [79] -21.862895 -21.976902 -22.091554 -22.193897 -22.277566 -22.340809
## [85] -22.382542 -22.400534 -22.393595 -22.365753 -22.328219 -22.296538
## [91] -22.284205Here we can see that the plots for hte realizations on death,positive and negative increase over time. This means that our variables here are cumulative to date and do not reflect the days’s total. We will need to focus on the day’s total as it will be difficult to forecast the cumulative increase.
# This plot shows the increase in positive cases over time.
plotts.sample.wge(us2$positiveIncrease)
## $autplt
## [1] 1.0000000 0.9699791 0.9416050 0.9101439 0.8855958 0.8589376 0.8400759
## [8] 0.8169588 0.7822600 0.7449305 0.7068916 0.6753458 0.6479670 0.6284668
## [15] 0.6008752 0.5682360 0.5368587 0.5059854 0.4751649 0.4501796 0.4270376
## [22] 0.4019983 0.3694803 0.3389772 0.3094579 0.2845992
##
## $freq
## [1] 0.005494505 0.010989011 0.016483516 0.021978022 0.027472527 0.032967033
## [7] 0.038461538 0.043956044 0.049450549 0.054945055 0.060439560 0.065934066
## [13] 0.071428571 0.076923077 0.082417582 0.087912088 0.093406593 0.098901099
## [19] 0.104395604 0.109890110 0.115384615 0.120879121 0.126373626 0.131868132
## [25] 0.137362637 0.142857143 0.148351648 0.153846154 0.159340659 0.164835165
## [31] 0.170329670 0.175824176 0.181318681 0.186813187 0.192307692 0.197802198
## [37] 0.203296703 0.208791209 0.214285714 0.219780220 0.225274725 0.230769231
## [43] 0.236263736 0.241758242 0.247252747 0.252747253 0.258241758 0.263736264
## [49] 0.269230769 0.274725275 0.280219780 0.285714286 0.291208791 0.296703297
## [55] 0.302197802 0.307692308 0.313186813 0.318681319 0.324175824 0.329670330
## [61] 0.335164835 0.340659341 0.346153846 0.351648352 0.357142857 0.362637363
## [67] 0.368131868 0.373626374 0.379120879 0.384615385 0.390109890 0.395604396
## [73] 0.401098901 0.406593407 0.412087912 0.417582418 0.423076923 0.428571429
## [79] 0.434065934 0.439560440 0.445054945 0.450549451 0.456043956 0.461538462
## [85] 0.467032967 0.472527473 0.478021978 0.483516484 0.489010989 0.494505495
## [91] 0.500000000
##
## $db
## [1] 14.2227365 15.5032402 7.5946898 9.3155600 4.3764892 -0.2424542
## [7] 2.7508484 -0.4922528 0.7737889 -1.8367628 -2.0630132 -3.2634822
## [13] -5.1059167 -7.1128370 -5.8394093 -7.4928164 -7.4106317 -7.4796411
## [19] -5.9926920 -7.9286420 -6.7417858 -7.3580975 -7.7715833 -11.4592878
## [25] -5.7130098 -4.0608460 -7.7314091 -17.3667018 -9.0950986 -10.5934114
## [31] -16.5277840 -10.7046757 -15.5827720 -12.4713924 -16.6740868 -13.7212208
## [37] -18.2049079 -15.0225542 -15.0033604 -13.0218576 -13.4832974 -10.9151165
## [43] -17.3205764 -11.1989911 -13.4534507 -13.5889044 -13.9459264 -15.8941852
## [49] -15.4243819 -20.0339780 -20.4442487 -10.6156562 -12.4719826 -10.7237344
## [55] -12.1267225 -18.6693407 -12.9421949 -17.6106540 -15.1012894 -23.0528491
## [61] -16.3859740 -19.8226319 -19.4515409 -19.1145075 -17.6059607 -21.0681024
## [67] -13.8648211 -17.8686977 -15.6460795 -21.3868724 -17.2343263 -18.3511836
## [73] -17.8696918 -20.9601880 -20.3303018 -18.9897721 -16.0003732 -14.9019851
## [79] -16.5756989 -14.1630850 -13.4348751 -16.3095403 -13.4669866 -14.1307274
## [85] -15.9978261 -12.1004576 -16.5118846 -12.2635868 -15.7331035 -19.3062757
## [91] -19.7004995
##
## $dbz
## [1] 11.9622903 11.6075420 11.0141762 10.1795761 9.1014793 7.7801350
## [7] 6.2224364 4.4493925 2.5081586 0.4866106 -1.4820724 -3.2496044
## [13] -4.7193340 -5.8832024 -6.7754329 -7.4094224 -7.7816484 -7.9218671
## [19] -7.9105581 -7.8425995 -7.7866132 -7.7743090 -7.8132538 -7.9056982
## [25] -8.0617099 -8.3019852 -8.6521629 -9.1339184 -9.7571082 -10.5143600
## [31] -11.3780099 -12.2998808 -13.2157071 -14.0555865 -14.7581785 -15.2826602
## [37] -15.6146048 -15.7672810 -15.7792516 -15.7055889 -15.6032574 -15.5169920
## [43] -15.4712881 -15.4692869 -15.4966339 -15.5285344 -15.5387810 -15.5090362
## [49] -15.4356889 -15.3318170 -15.2238891 -15.1452705 -15.1293821 -15.2043203
## [55] -15.3892464 -15.6919481 -16.1068741 -16.6134394 -17.1753881 -17.7432574
## [61] -18.2624631 -18.6871772 -18.9946082 -19.1908666 -19.3045915 -19.3735425
## [67] -19.4319508 -19.5020043 -19.5887229 -19.6772576 -19.7339371 -19.7143633
## [73] -19.5799089 -19.3162504 -18.9412401 -18.4960382 -18.0277404 -17.5758973
## [79] -17.1677268 -16.8195809 -16.5406040 -16.3358975 -16.2081069 -16.1573605
## [85] -16.1799827 -16.2665968 -16.4003572 -16.5562924 -16.7030465 -16.8081434
## [91] -16.8463737# This shows the same for negative cases over time.
plotts.sample.wge(us2$negativeIncrease)
## $autplt
## [1] 1.0000000 0.9688031 0.9515967 0.9329244 0.9173372 0.9008759 0.8851659
## [8] 0.8717487 0.8518122 0.8331232 0.8136295 0.8004597 0.7849694 0.7749357
## [15] 0.7643802 0.7480473 0.7329418 0.7157441 0.7027817 0.6857972 0.6709983
## [22] 0.6568792 0.6388638 0.6189034 0.6001762 0.5859067
##
## $freq
## [1] 0.005494505 0.010989011 0.016483516 0.021978022 0.027472527 0.032967033
## [7] 0.038461538 0.043956044 0.049450549 0.054945055 0.060439560 0.065934066
## [13] 0.071428571 0.076923077 0.082417582 0.087912088 0.093406593 0.098901099
## [19] 0.104395604 0.109890110 0.115384615 0.120879121 0.126373626 0.131868132
## [25] 0.137362637 0.142857143 0.148351648 0.153846154 0.159340659 0.164835165
## [31] 0.170329670 0.175824176 0.181318681 0.186813187 0.192307692 0.197802198
## [37] 0.203296703 0.208791209 0.214285714 0.219780220 0.225274725 0.230769231
## [43] 0.236263736 0.241758242 0.247252747 0.252747253 0.258241758 0.263736264
## [49] 0.269230769 0.274725275 0.280219780 0.285714286 0.291208791 0.296703297
## [55] 0.302197802 0.307692308 0.313186813 0.318681319 0.324175824 0.329670330
## [61] 0.335164835 0.340659341 0.346153846 0.351648352 0.357142857 0.362637363
## [67] 0.368131868 0.373626374 0.379120879 0.384615385 0.390109890 0.395604396
## [73] 0.401098901 0.406593407 0.412087912 0.417582418 0.423076923 0.428571429
## [79] 0.434065934 0.439560440 0.445054945 0.450549451 0.456043956 0.461538462
## [85] 0.467032967 0.472527473 0.478021978 0.483516484 0.489010989 0.494505495
## [91] 0.500000000
##
## $db
## [1] 17.8786349 9.7933404 6.6915875 5.4268246 2.9044749 0.5613665
## [7] 0.3314022 -2.0749363 0.1948281 -3.0823919 -2.3477911 -3.7581763
## [13] -4.1464510 -10.2266176 -8.0028120 -7.6605968 -9.6781694 -9.7372330
## [19] -6.2172693 -11.4200424 -6.9059420 -13.2480454 -11.7793943 -8.5162704
## [25] -5.8584139 -9.9622266 -14.6340579 -10.8485900 -8.6456673 -15.6863451
## [31] -13.4759449 -12.3040982 -14.1727577 -10.2813822 -12.3303527 -12.1203098
## [37] -16.3776747 -15.7282307 -11.9044288 -18.1352683 -12.6460095 -10.9069671
## [43] -17.7261080 -12.7269645 -13.5853925 -21.9701066 -14.3325474 -12.3474954
## [49] -10.7210470 -21.1158795 -11.6743587 -10.4159770 -14.7304017 -16.3735978
## [55] -15.6092644 -23.3224542 -12.2152153 -29.6853921 -15.4559060 -19.5583003
## [61] -12.9335015 -27.1092661 -24.2818225 -13.3118464 -13.1328630 -18.2829810
## [67] -13.5133307 -29.1606958 -11.6244503 -18.4211246 -13.1254641 -29.1177110
## [73] -19.7062042 -16.4823777 -23.6088298 -19.9097553 -13.6295672 -17.8423222
## [79] -10.0594879 -14.9020687 -12.6025210 -20.7907400 -11.3684985 -14.2020457
## [85] -14.2903242 -20.1120872 -18.4127696 -26.2823412 -12.7005081 -15.4224356
## [91] -19.6614861
##
## $dbz
## [1] 12.2411091 11.8452204 11.1808831 10.2417877 9.0204888 7.5111146
## [7] 5.7156510 3.6576577 1.4086663 -0.8758283 -2.9548502 -4.6210632
## [13] -5.8722986 -6.8672503 -7.7151184 -8.3853322 -8.7942916 -8.9485615
## [19] -8.9742836 -9.0187503 -9.1659373 -9.4226071 -9.7433056 -10.0692974
## [25] -10.3675149 -10.6485597 -10.9518528 -11.3135536 -11.7418266 -12.2105637
## [31] -12.6716636 -13.0805112 -13.4199128 -13.7037136 -13.9578980 -14.1976442
## [37] -14.4184064 -14.6043153 -14.7444825 -14.8429827 -14.9144888 -14.9707913
## [43] -15.0105443 -15.0204601 -14.9871451 -14.9109153 -14.8110882 -14.7194383
## [49] -14.6680996 -14.6801969 -14.7665262 -14.9268559 -15.1529187 -15.4307215
## [55] -15.7411281 -16.0593658 -16.3556913 -16.5997660 -16.7689476 -16.8563946
## [61] -16.8729551 -16.8409994 -16.7848307 -16.7239361 -16.6715762 -16.6369103
## [67] -16.6271478 -16.6471122 -16.6958311 -16.7619620 -16.8212753 -16.8397786
## [73] -16.7843310 -16.6374546 -16.4073653 -16.1257843 -15.8357556 -15.5784621
## [79] -15.3853916 -15.2761602 -15.2593359 -15.3338597 -15.4899638 -15.7096535
## [85] -15.9676294 -16.2338804 -16.4786298 -16.6785621 -16.8212311 -16.9045897
## [91] -16.9317276The totalIncrease created by summing the positive and negative increase should be equivalent to the number of tests administered without considering pending cases.
totalIncrease = us2$positiveIncrease + us2$negativeIncrease
# This variable (totalTestResultsIncresase) should closely match the totalIncrease
plotts.sample.wge(totalIncrease)
## $autplt
## [1] 1.0000000 0.9697066 0.9515991 0.9319392 0.9154431 0.8983689 0.8821088
## [8] 0.8681424 0.8471479 0.8275953 0.8069882 0.7931058 0.7771362 0.7665773
## [15] 0.7550266 0.7380453 0.7226435 0.7047302 0.6910972 0.6740446 0.6591417
## [22] 0.6445669 0.6261366 0.6055607 0.5867574 0.5722355
##
## $freq
## [1] 0.005494505 0.010989011 0.016483516 0.021978022 0.027472527 0.032967033
## [7] 0.038461538 0.043956044 0.049450549 0.054945055 0.060439560 0.065934066
## [13] 0.071428571 0.076923077 0.082417582 0.087912088 0.093406593 0.098901099
## [19] 0.104395604 0.109890110 0.115384615 0.120879121 0.126373626 0.131868132
## [25] 0.137362637 0.142857143 0.148351648 0.153846154 0.159340659 0.164835165
## [31] 0.170329670 0.175824176 0.181318681 0.186813187 0.192307692 0.197802198
## [37] 0.203296703 0.208791209 0.214285714 0.219780220 0.225274725 0.230769231
## [43] 0.236263736 0.241758242 0.247252747 0.252747253 0.258241758 0.263736264
## [49] 0.269230769 0.274725275 0.280219780 0.285714286 0.291208791 0.296703297
## [55] 0.302197802 0.307692308 0.313186813 0.318681319 0.324175824 0.329670330
## [61] 0.335164835 0.340659341 0.346153846 0.351648352 0.357142857 0.362637363
## [67] 0.368131868 0.373626374 0.379120879 0.384615385 0.390109890 0.395604396
## [73] 0.401098901 0.406593407 0.412087912 0.417582418 0.423076923 0.428571429
## [79] 0.434065934 0.439560440 0.445054945 0.450549451 0.456043956 0.461538462
## [85] 0.467032967 0.472527473 0.478021978 0.483516484 0.489010989 0.494505495
## [91] 0.500000000
##
## $db
## [1] 17.7282550 10.3240493 6.8437548 5.8713444 3.1103259 0.5689188
## [7] 0.6194258 -1.8679340 0.3214384 -2.8973864 -2.2408720 -3.6361499
## [13] -4.1303870 -9.8715017 -7.7414080 -7.5614942 -9.3998914 -9.4727763
## [19] -6.1136728 -11.0701121 -6.8136620 -12.5893425 -11.4105647 -8.9573854
## [25] -5.7723114 -9.3327707 -13.8961501 -11.1236411 -8.6370668 -15.1014398
## [31] -13.5839541 -12.0969930 -14.1853223 -10.3515792 -12.5122408 -12.1985407
## [37] -16.4793229 -15.5887339 -12.3803328 -19.0592378 -12.6437405 -10.8224995
## [43] -17.7757809 -12.5809730 -13.5756202 -21.0178315 -14.3013267 -12.5237212
## [49] -10.9301170 -20.9433502 -12.1168100 -10.4945093 -15.0894533 -15.7113619
## [55] -15.3516567 -23.7518674 -12.1840463 -29.2847164 -15.3424045 -20.2736452
## [61] -13.2490245 -26.8699947 -23.7421991 -13.7469668 -13.4218223 -18.3778469
## [67] -13.5273675 -27.5621762 -11.7853893 -18.5354788 -13.4555084 -31.9611663
## [73] -19.9597872 -16.9292636 -23.2794616 -19.9126294 -13.7034345 -17.5543916
## [79] -10.3228936 -14.8542849 -12.5936968 -20.2925758 -11.4290997 -14.1541159
## [85] -14.5691480 -20.6246748 -18.3072661 -26.1116726 -13.1523594 -15.5738143
## [91] -19.5776324
##
## $dbz
## [1] 12.2257976 11.8321166 11.1716838 10.2385665 9.0259278 7.5288155
## [7] 5.7504929 3.7159584 1.4972824 -0.7532946 -2.8049002 -4.4613191
## [13] -5.7184635 -6.7228544 -7.5766884 -8.2516569 -8.6693329 -8.8363964
## [19] -8.8736573 -8.9234765 -9.0682002 -9.3157310 -9.6236993 -9.9375395
## [25] -10.2277992 -10.5066105 -10.8130806 -11.1828920 -11.6249140 -12.1149084
## [31] -12.6064144 -13.0542963 -13.4380131 -13.7667996 -14.0628318 -14.3387243
## [37] -14.5877453 -14.7917564 -14.9383355 -15.0320456 -15.0904360 -15.1297917
## [43] -15.1537736 -15.1535492 -15.1184645 -15.0489196 -14.9617539 -14.8849335
## [49] -14.8468925 -14.8680084 -14.9575462 -15.1150231 -15.3332653 -15.6007208
## [55] -15.9018144 -16.2158468 -16.5165955 -16.7752913 -16.9677992 -17.0828853
## [61] -17.1258616 -17.1146779 -17.0716352 -17.0166097 -16.9648667 -16.9282482
## [67] -16.9164697 -16.9360006 -16.9861995 -17.0545423 -17.1142453 -17.1280916
## [73] -17.0604623 -16.8937046 -16.6385762 -16.3307943 -16.0170673 -15.7409801
## [79] -15.5352554 -15.4200956 -15.4044539 -15.4876685 -15.6603520 -15.9046158
## [85] -16.1945359 -16.4982206 -16.7824912 -17.0194922 -17.1921545 -17.2948538
## [91] -17.3286462The totalIncrease variable shows to match exactly with the totalTestResultsIncrease variable, so now we can have confidence on how that variable was calculated and use it moving forward.
# Here is a graph of increase count in total tests.
plotts.sample.wge(us2$totalTestResultsIncrease)
## $autplt
## [1] 1.0000000 0.9697066 0.9515991 0.9319392 0.9154431 0.8983689 0.8821088
## [8] 0.8681424 0.8471479 0.8275953 0.8069882 0.7931058 0.7771362 0.7665773
## [15] 0.7550266 0.7380453 0.7226435 0.7047302 0.6910972 0.6740446 0.6591417
## [22] 0.6445669 0.6261366 0.6055607 0.5867574 0.5722355
##
## $freq
## [1] 0.005494505 0.010989011 0.016483516 0.021978022 0.027472527 0.032967033
## [7] 0.038461538 0.043956044 0.049450549 0.054945055 0.060439560 0.065934066
## [13] 0.071428571 0.076923077 0.082417582 0.087912088 0.093406593 0.098901099
## [19] 0.104395604 0.109890110 0.115384615 0.120879121 0.126373626 0.131868132
## [25] 0.137362637 0.142857143 0.148351648 0.153846154 0.159340659 0.164835165
## [31] 0.170329670 0.175824176 0.181318681 0.186813187 0.192307692 0.197802198
## [37] 0.203296703 0.208791209 0.214285714 0.219780220 0.225274725 0.230769231
## [43] 0.236263736 0.241758242 0.247252747 0.252747253 0.258241758 0.263736264
## [49] 0.269230769 0.274725275 0.280219780 0.285714286 0.291208791 0.296703297
## [55] 0.302197802 0.307692308 0.313186813 0.318681319 0.324175824 0.329670330
## [61] 0.335164835 0.340659341 0.346153846 0.351648352 0.357142857 0.362637363
## [67] 0.368131868 0.373626374 0.379120879 0.384615385 0.390109890 0.395604396
## [73] 0.401098901 0.406593407 0.412087912 0.417582418 0.423076923 0.428571429
## [79] 0.434065934 0.439560440 0.445054945 0.450549451 0.456043956 0.461538462
## [85] 0.467032967 0.472527473 0.478021978 0.483516484 0.489010989 0.494505495
## [91] 0.500000000
##
## $db
## [1] 17.7282550 10.3240493 6.8437548 5.8713444 3.1103259 0.5689188
## [7] 0.6194258 -1.8679340 0.3214384 -2.8973864 -2.2408720 -3.6361499
## [13] -4.1303870 -9.8715017 -7.7414080 -7.5614942 -9.3998914 -9.4727763
## [19] -6.1136728 -11.0701121 -6.8136620 -12.5893425 -11.4105647 -8.9573854
## [25] -5.7723114 -9.3327707 -13.8961501 -11.1236411 -8.6370668 -15.1014398
## [31] -13.5839541 -12.0969930 -14.1853223 -10.3515792 -12.5122408 -12.1985407
## [37] -16.4793229 -15.5887339 -12.3803328 -19.0592378 -12.6437405 -10.8224995
## [43] -17.7757809 -12.5809730 -13.5756202 -21.0178315 -14.3013267 -12.5237212
## [49] -10.9301170 -20.9433502 -12.1168100 -10.4945093 -15.0894533 -15.7113619
## [55] -15.3516567 -23.7518674 -12.1840463 -29.2847164 -15.3424045 -20.2736452
## [61] -13.2490245 -26.8699947 -23.7421991 -13.7469668 -13.4218223 -18.3778469
## [67] -13.5273675 -27.5621762 -11.7853893 -18.5354788 -13.4555084 -31.9611663
## [73] -19.9597872 -16.9292636 -23.2794616 -19.9126294 -13.7034345 -17.5543916
## [79] -10.3228936 -14.8542849 -12.5936968 -20.2925758 -11.4290997 -14.1541159
## [85] -14.5691480 -20.6246748 -18.3072661 -26.1116726 -13.1523594 -15.5738143
## [91] -19.5776324
##
## $dbz
## [1] 12.2257976 11.8321166 11.1716838 10.2385665 9.0259278 7.5288155
## [7] 5.7504929 3.7159584 1.4972824 -0.7532946 -2.8049002 -4.4613191
## [13] -5.7184635 -6.7228544 -7.5766884 -8.2516569 -8.6693329 -8.8363964
## [19] -8.8736573 -8.9234765 -9.0682002 -9.3157310 -9.6236993 -9.9375395
## [25] -10.2277992 -10.5066105 -10.8130806 -11.1828920 -11.6249140 -12.1149084
## [31] -12.6064144 -13.0542963 -13.4380131 -13.7667996 -14.0628318 -14.3387243
## [37] -14.5877453 -14.7917564 -14.9383355 -15.0320456 -15.0904360 -15.1297917
## [43] -15.1537736 -15.1535492 -15.1184645 -15.0489196 -14.9617539 -14.8849335
## [49] -14.8468925 -14.8680084 -14.9575462 -15.1150231 -15.3332653 -15.6007208
## [55] -15.9018144 -16.2158468 -16.5165955 -16.7752913 -16.9677992 -17.0828853
## [61] -17.1258616 -17.1146779 -17.0716352 -17.0166097 -16.9648667 -16.9282482
## [67] -16.9164697 -16.9360006 -16.9861995 -17.0545423 -17.1142453 -17.1280916
## [73] -17.0604623 -16.8937046 -16.6385762 -16.3307943 -16.0170673 -15.7409801
## [79] -15.5352554 -15.4200956 -15.4044539 -15.4876685 -15.6603520 -15.9046158
## [85] -16.1945359 -16.4982206 -16.7824912 -17.0194922 -17.1921545 -17.2948538
## [91] -17.3286462Now we are divding the total test by cumulative positive cases.
pp = totalIncrease / us2$positive
plotts.sample.wge(pp)
## $autplt
## [1] 1.000000000 0.869707176 0.840284240 0.830714964 0.786400348
## [6] 0.769950222 0.725927290 0.702261776 0.659974660 0.617738778
## [11] 0.518488665 0.469354845 0.456471593 0.382306058 0.354189790
## [16] 0.287449638 0.236542805 0.190566212 0.134284329 0.096891654
## [21] 0.041937922 0.015952649 -0.002207235 -0.042464952 -0.075364873
## [26] -0.100435187
##
## $freq
## [1] 0.005494505 0.010989011 0.016483516 0.021978022 0.027472527 0.032967033
## [7] 0.038461538 0.043956044 0.049450549 0.054945055 0.060439560 0.065934066
## [13] 0.071428571 0.076923077 0.082417582 0.087912088 0.093406593 0.098901099
## [19] 0.104395604 0.109890110 0.115384615 0.120879121 0.126373626 0.131868132
## [25] 0.137362637 0.142857143 0.148351648 0.153846154 0.159340659 0.164835165
## [31] 0.170329670 0.175824176 0.181318681 0.186813187 0.192307692 0.197802198
## [37] 0.203296703 0.208791209 0.214285714 0.219780220 0.225274725 0.230769231
## [43] 0.236263736 0.241758242 0.247252747 0.252747253 0.258241758 0.263736264
## [49] 0.269230769 0.274725275 0.280219780 0.285714286 0.291208791 0.296703297
## [55] 0.302197802 0.307692308 0.313186813 0.318681319 0.324175824 0.329670330
## [61] 0.335164835 0.340659341 0.346153846 0.351648352 0.357142857 0.362637363
## [67] 0.368131868 0.373626374 0.379120879 0.384615385 0.390109890 0.395604396
## [73] 0.401098901 0.406593407 0.412087912 0.417582418 0.423076923 0.428571429
## [79] 0.434065934 0.439560440 0.445054945 0.450549451 0.456043956 0.461538462
## [85] 0.467032967 0.472527473 0.478021978 0.483516484 0.489010989 0.494505495
## [91] 0.500000000
##
## $db
## [1] 14.126881 14.281246 11.561975 8.018969 6.099023 2.516908
## [7] -4.496788 -6.398092 -5.504191 -21.819520 -6.791808 -5.762722
## [13] -8.821942 -10.343204 -12.878167 -8.152340 -8.906452 -11.043667
## [19] -10.843025 -14.110123 -7.592559 -6.131387 -4.999954 -3.847412
## [25] -3.300102 -6.095248 -9.628748 -8.955256 -5.962344 -7.920908
## [31] -11.278160 -13.031733 -15.216341 -17.811117 -20.117617 -23.630139
## [37] -13.779093 -11.442714 -13.530758 -8.351224 -4.729573 -3.455995
## [43] -9.228062 -16.444414 -7.272109 -5.467860 -5.776384 -11.388356
## [49] -9.563139 -7.670562 -11.715271 -13.583353 -17.937904 -13.953516
## [55] -12.322573 -15.460384 -12.229460 -5.827924 -5.309766 -8.228840
## [61] -9.838321 -5.971913 -5.611708 -6.761064 -9.632216 -6.673955
## [67] -6.153850 -12.124588 -20.626310 -20.163589 -28.144717 -19.286203
## [73] -11.341492 -5.875255 -6.475635 -7.397755 -5.410746 -5.666943
## [79] -7.591386 -8.847660 -10.104192 -9.412369 -13.244141 -17.190576
## [85] -20.455397 -14.272991 -17.520270 -9.819296 -10.097260 -16.908415
## [91] -30.071673
##
## $dbz
## [1] 11.2681884 10.9787369 10.4933831 9.8077717 8.9159696 7.8108389
## [7] 6.4849633 4.9328157 3.1556558 1.1721784 -0.9600851 -3.1098856
## [13] -5.0458937 -6.4976211 -7.3374274 -7.6661923 -7.6682794 -7.4770427
## [19] -7.1721138 -6.8165855 -6.4711285 -6.1884317 -6.0063172 -5.9469059
## [25] -6.0201787 -6.2286456 -6.5710726 -7.0442472 -7.6420135 -8.3505086
## [31] -9.1380820 -9.9392537 -10.6390412 -11.0825447 -11.1428436 -10.8142150
## [37] -10.2192217 -9.5194688 -8.8403639 -8.2553614 -7.7995933 -7.4857260
## [43] -7.3145153 -7.2807189 -7.3761535 -7.5909825 -7.9135103 -8.3281717
## [49] -8.8111595 -9.3236548 -9.8049392 -10.1727128 -10.3420045 -10.2636639
## [55] -9.9552310 -9.4918803 -8.9671196 -8.4597960 -8.0231645 -7.6883601
## [61] -7.4717623 -7.3811069 -7.4189837 -7.5836218 -7.8669394 -8.2496824
## [67] -8.6938975 -9.1353872 -9.4845887 -9.6492094 -9.5791162 -9.2999418
## [73] -8.8989556 -8.4791915 -8.1245656 -7.8901485 -7.8066710 -7.8877740
## [79] -8.1351896 -8.5409651 -9.0870265 -9.7426770 -10.4613142 -11.1794630
## [85] -11.8235574 -12.3286447 -12.6629989 -12.8404735 -12.9085897 -12.9223799
## [91] -12.9221699plotts.sample.wge(pp[40:140])
## $autplt
## [1] 1.000000000 0.849145416 0.817853490 0.809471859 0.771557744
## [6] 0.763231902 0.720605060 0.699006876 0.637690444 0.598835264
## [11] 0.492085673 0.447179938 0.446140305 0.376994375 0.346652286
## [16] 0.278435422 0.233806187 0.199898660 0.153951516 0.134610290
## [21] 0.092668556 0.075405189 0.067186162 0.024926953 0.003291947
## [26] -0.024498240
##
## $freq
## [1] 0.00990099 0.01980198 0.02970297 0.03960396 0.04950495 0.05940594
## [7] 0.06930693 0.07920792 0.08910891 0.09900990 0.10891089 0.11881188
## [13] 0.12871287 0.13861386 0.14851485 0.15841584 0.16831683 0.17821782
## [19] 0.18811881 0.19801980 0.20792079 0.21782178 0.22772277 0.23762376
## [25] 0.24752475 0.25742574 0.26732673 0.27722772 0.28712871 0.29702970
## [31] 0.30693069 0.31683168 0.32673267 0.33663366 0.34653465 0.35643564
## [37] 0.36633663 0.37623762 0.38613861 0.39603960 0.40594059 0.41584158
## [43] 0.42574257 0.43564356 0.44554455 0.45544554 0.46534653 0.47524752
## [49] 0.48514851 0.49504950
##
## $db
## [1] 14.452829 10.762374 5.119802 -2.526528 -7.288997 -14.102895
## [7] -18.310311 -9.324752 -7.440743 -8.210643 -18.989978 -6.941506
## [13] -4.420881 -2.270410 -6.450405 -6.936395 -8.823309 -12.669066
## [19] -16.437682 -18.341215 -15.380102 -6.437235 -6.251663 -8.046733
## [25] -6.480863 -6.992158 -6.917435 -10.210531 -13.193867 -14.506619
## [31] -11.066232 -7.345723 -4.854577 -7.369940 -5.644268 -6.057265
## [37] -7.813122 -13.035399 -30.163794 -17.475631 -6.527271 -4.959161
## [43] -4.394756 -8.641767 -9.792611 -9.209634 -14.550936 -23.856124
## [49] -11.503141 -16.479205
##
## $dbz
## [1] 10.2820019 9.6282061 8.5264169 6.9589589 4.9038928 2.3452826
## [7] -0.6872321 -3.9761143 -6.7699098 -7.9322706 -7.5876059 -6.7926032
## [13] -6.1205069 -5.7565048 -5.7566186 -6.1175469 -6.7749203 -7.5739273
## [19] -8.2506896 -8.5290644 -8.3498423 -7.9252859 -7.5241517 -7.3160007
## [25] -7.3629819 -7.6469128 -8.0784187 -8.4986964 -8.7197809 -8.6333122
## [31] -8.2970713 -7.8817610 -7.5544942 -7.4201929 -7.5157511 -7.8104350
## [37] -8.1966031 -8.4963520 -8.5409805 -8.3060438 -7.9388680 -7.6379165
## [43] -7.5488669 -7.7453025 -8.2438629 -9.0141140 -9.9745991 -10.9825850
## [49] -11.8389549 -12.3350874Now we are creating the variable posper, which stands for posive percentage, given as decimal percentages.Since our positive percentage varies widely with little data in the first 40 days of the dataset, so we will focus on the data from day 40 to 140.
posper = (us2$positiveIncrease[50:140]/us2$totalTestResultsIncrease[50:140])
#posper = (us2$positiveIncrease/us2$totalTestResultsIncrease)
plotts.sample.wge(posper)
## $autplt
## [1] 1.0000000 0.8708698 0.8746167 0.8704957 0.8210217 0.8290970 0.7950567
## [8] 0.7616090 0.7566081 0.7297009 0.6975731 0.6484219 0.6352532 0.5861958
## [15] 0.5688315 0.5338068 0.4769021 0.4559543 0.4085230 0.3872615 0.3238917
## [22] 0.2732022 0.2404804 0.1995276 0.1478254 0.1233541
##
## $freq
## [1] 0.01098901 0.02197802 0.03296703 0.04395604 0.05494505 0.06593407
## [7] 0.07692308 0.08791209 0.09890110 0.10989011 0.12087912 0.13186813
## [13] 0.14285714 0.15384615 0.16483516 0.17582418 0.18681319 0.19780220
## [19] 0.20879121 0.21978022 0.23076923 0.24175824 0.25274725 0.26373626
## [25] 0.27472527 0.28571429 0.29670330 0.30769231 0.31868132 0.32967033
## [31] 0.34065934 0.35164835 0.36263736 0.37362637 0.38461538 0.39560440
## [37] 0.40659341 0.41758242 0.42857143 0.43956044 0.45054945 0.46153846
## [43] 0.47252747 0.48351648 0.49450549
##
## $db
## [1] 15.857855 -1.437845 -4.196407 -4.975054 -3.842813 -3.261753
## [7] -16.021221 -17.334291 -8.247654 -3.378722 -9.879409 -31.943441
## [13] -6.599601 -22.042356 -17.934454 -15.304123 -11.882936 -10.189475
## [19] -9.749542 -7.261349 -13.966675 -13.914466 -13.282440 -13.544477
## [25] -8.199119 -13.837319 -14.961802 -18.695208 -4.703262 -12.493558
## [31] -5.687354 -3.313666 -9.517705 -5.286433 -14.803735 -7.465879
## [37] -4.620027 -10.330086 -8.687790 -13.048557 -20.323370 -15.235530
## [43] -15.277547 -4.855966 -20.358868
##
## $dbz
## [1] 10.4352044 9.6030940 8.1940306 6.1784175 3.5328617 0.3025509
## [7] -3.1985694 -6.0963581 -7.6114782 -8.2090983 -8.5282311 -8.7520179
## [13] -8.9881595 -9.3573513 -9.8399907 -10.2254314 -10.2983043 -10.0973420
## [19] -9.8558105 -9.7592977 -9.8654420 -10.1459305 -10.5263399 -10.8865816
## [25] -11.0534973 -10.8648459 -10.3045847 -9.5313643 -8.7480211 -8.0926493
## [31] -7.6228900 -7.3382892 -7.2043527 -7.1765019 -7.2263848 -7.3605123
## [37] -7.6160099 -8.0326046 -8.6143299 -9.2940737 -9.9215782 -10.3218861
## [43] -10.4288439 -10.3473416 -10.2503127The realization of the posper plot shows an decreasing trend over time for approximately the past sixty days of data indicating that positive test results are decreasing relative to the number of administered tests.
Proceed with stationary model being aware that short term stationarity is possibly unlikely.
aic5.wge(posper, p=0:5, q=0:2) # arma(2,2)## ---------WORKING... PLEASE WAIT...
##
##
## Five Smallest Values of aic## p q aic
## 9 2 2 -7.450198
## 10 3 0 -7.415703
## 16 5 0 -7.394922
## 11 3 1 -7.394331
## 13 4 0 -7.393858aic5.wge(posper, p=0:5, q=0:2, type='bic') # arma(2,2)## ---------WORKING... PLEASE WAIT...
##
##
## Five Smallest Values of bic## p q bic
## 9 2 2 -7.312238
## 10 3 0 -7.305335
## 5 1 1 -7.281966
## 8 2 1 -7.274844
## 6 1 2 -7.265863Estimate the model for the p = 2 and q = 2
estPosper22 = est.arma.wge(posper, p=2, q=2)##
## Coefficients of Original polynomial:
## 1.9642 -0.9677
##
## Factor Roots Abs Recip System Freq
## 1-1.9642B+0.9677B^2 1.0149+-0.0575i 0.9837 0.0090
##
## mean(posper)## [1] 0.125892estPosper22$phi## [1] 1.9642393 -0.9676536estPosper22$theta## [1] 1.8489656 -0.8875016Based on our estimation, the final model for positive percentage is:
(1 - 1.96B + 0.97B^2)(X - 0.13) = (1 -1.185a_t + 0.89a_t^2)
Create the Factor tables for the estimated model we created.
factor.wge(phi=estPosper22$phi)##
## Coefficients of Original polynomial:
## 1.9642 -0.9677
##
## Factor Roots Abs Recip System Freq
## 1-1.9642B+0.9677B^2 1.0149+-0.0575i 0.9837 0.0090
##
## Forecasting
Forcasting 10 ahead for both ARMA models. Both of these will eventually reach the mean of 0.16
fore.arma.wge(posper,phi=estPosper22$phi, theta=estPosper22$theta, n.ahead=10,plot=TRUE)
## $f
## [1] 0.05170743 0.05177394 0.05209136 0.05265049 0.05344160 0.05445449
## [7] 0.05567853 0.05710270 0.05871568 0.06050584
##
## $ll
## [1] -0.001004122 -0.001286676 -0.001526527 -0.001762016 -0.002020339
## [6] -0.002317712 -0.002660313 -0.003045834 -0.003465404 -0.003905594
##
## $ul
## [1] 0.1044190 0.1048346 0.1057092 0.1070630 0.1089035 0.1112267 0.1140174
## [8] 0.1172512 0.1208968 0.1249173
##
## $resid
## [1] 0.0000000000 0.0000000000 -0.0073797183 0.0891554359 0.0537685407
## [6] -0.0236956311 0.1038957821 0.0141749394 0.0532589863 0.0517683006
## [11] 0.0257323541 0.0715921667 0.0528920948 -0.0018723369 -0.0099540774
## [16] 0.0054326860 0.0109388846 0.0029467473 0.0366425047 -0.0125181161
## [21] 0.0177424318 0.0280375193 0.0174160054 0.0144652063 -0.0883315987
## [26] -0.0154257662 -0.0392272027 -0.0272901829 -0.0162690676 -0.0188969957
## [31] 0.0010420787 0.0100777453 -0.0106384335 -0.0195668932 -0.0322452513
## [36] 0.0280256738 -0.0022455478 0.0151814170 0.0083990256 -0.0016467927
## [41] 0.0016910809 0.0003098111 -0.0785621464 0.0100250867 -0.0045761866
## [46] -0.0128731426 -0.0025990574 -0.0130321785 0.0010385872 -0.0033875733
## [51] 0.0190622719 0.0074793141 0.0178611777 0.0112370946 0.0021085539
## [56] -0.0095295324 0.0141499919 0.0043249031 0.0077095185 0.0012047599
## [61] -0.0017364091 -0.0332555948 -0.0025372783 -0.0072714440 0.0023219264
## [66] -0.0005224882 0.0002792402 -0.0123586636 -0.0050923385 -0.0105738475
## [71] -0.0078937665 -0.0008013121 0.0030179131 -0.0011277722 -0.0031503357
## [76] -0.0107162401 0.0005265607 0.0091819180 -0.0027557348 -0.0067591250
## [81] 0.0008393554 -0.0011705865 -0.0071086394 -0.0092036972 -0.0134600402
## [86] -0.0112432954 -0.0098281878 -0.0103846957 -0.0127593901 -0.0109142816
## [91] -0.0119992790
##
## $wnv
## [1] 0.0007232685
##
## $se
## [1] 0.02689365 0.02707174 0.02735606 0.02776148 0.02829691 0.02896541
## [7] 0.02976471 0.03068803 0.03172504 0.03286298
##
## $psi
## [1] 0.1152737 0.1462732 0.1757706 0.2037137 0.2300573 0.2547634 0.2778005
## [8] 0.2991439 0.3187756 0.3366838Calculate the ASE for 10 steps back using the ARMA Models
for22 = fore.arma.wge(posper,phi=estPosper22$phi, theta=estPosper22$theta, n.ahead=10,lastn=TRUE, plot=TRUE)
ASE22 = mean((posper[(91-10+1):91] - for22$f)^2)
ASE22## [1] 0.0003449402For better confirmation of the performance of the model, we will run a windowed ASE which will test the performance of the model during 10 unit windows throughout the time series of the realization.
trainingSize = 70
horizon = 10
ASEHolder = numeric()
for( i in 1:(91-(trainingSize + horizon) + 1))
{
forecasts = fore.aruma.wge(posper[i:(i+(trainingSize-1))],phi = estPosper22$phi, theta = estPosper22$theta, s = 0, d = 0,n.ahead = horizon)
ASE = mean((posper[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
ASEHolder[i] = ASE
}
ASEHolder## [1] 9.436351e-05 4.880261e-05 3.108709e-04 4.371886e-05 3.146104e-05
## [6] 6.253707e-04 4.195351e-05 2.520207e-04 2.412527e-04 2.092395e-04
## [11] 5.091830e-04 3.328257e-04hist(ASEHolder)
WindowedASE = mean(ASEHolder)
summary(ASEHolder)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.146e-05 4.753e-05 2.252e-04 2.284e-04 3.164e-04 6.254e-04WindowedASE## [1] 0.0002284219Here, we can see that our models hover towards zero with our histogram suggesting our model is regularly hitting an ASE of 0.0009 and below.
Now we will look 20 steps back
for22 = fore.arma.wge(posper,phi=estPosper22$phi, theta=estPosper22$theta, n.ahead=20,lastn=TRUE, plot=TRUE)
ASE22 = mean((posper[(91-20+1):91] - for22$f)^2)
ASE22## [1] 0.0002598811Running the 20 day lookback forecast, we cans see that our model is predicting the positive percentage to increase back towards the mean which will ultimately oscillate if we extend the mean far out enough. Doing so on a larger window back, we ended up with a better ASE than the 10 day lookback window.
trainingSize = 70
horizon = 10
ASEHolder = numeric()
for( i in 1:(81-(trainingSize + horizon) + 1))
{
forecasts = fore.aruma.wge(posper[i:(i+(trainingSize-1))],phi = estPosper22$phi, theta = estPosper22$theta, s = 0, d = 0,n.ahead = horizon)
ASE = mean((posper[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
ASEHolder[i] = ASE
}
ASEHolder## [1] 9.436351e-05 4.880261e-05hist(ASEHolder)
WindowedASE = mean(ASEHolder)
summary(ASEHolder)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.880e-05 6.019e-05 7.158e-05 7.158e-05 8.297e-05 9.436e-05WindowedASE## [1] 7.158306e-05Now lets look into a seasonal model and see if there is a significant change in our ASE score.
Trend Modeling
p.ns = artrans.wge(posper,phi.tr=(1))
plotts.sample.wge(p.ns)
## $autplt
## [1] 1.00000000 -0.55695521 0.03373231 0.17963754 -0.23793813 0.18661135
## [7] -0.02340134 -0.10506081 0.08676206 0.02484071 0.07412824 -0.16828997
## [13] 0.15373153 -0.11793145 0.06468540 0.07756060 -0.14680375 0.11237394
## [19] -0.11345250 0.17835694 -0.06567108 -0.08112372 0.03954708 0.03713005
## [25] -0.08873895 0.09334638
##
## $freq
## [1] 0.01111111 0.02222222 0.03333333 0.04444444 0.05555556 0.06666667
## [7] 0.07777778 0.08888889 0.10000000 0.11111111 0.12222222 0.13333333
## [13] 0.14444444 0.15555556 0.16666667 0.17777778 0.18888889 0.20000000
## [19] 0.21111111 0.22222222 0.23333333 0.24444444 0.25555556 0.26666667
## [25] 0.27777778 0.28888889 0.30000000 0.31111111 0.32222222 0.33333333
## [31] 0.34444444 0.35555556 0.36666667 0.37777778 0.38888889 0.40000000
## [37] 0.41111111 0.42222222 0.43333333 0.44444444 0.45555556 0.46666667
## [43] 0.47777778 0.48888889 0.50000000
##
## $db
## [1] -4.9213122 -11.3931048 -23.3631652 -17.4268252 -14.9551672 -8.8762207
## [7] -14.2590846 -12.7082065 -8.5955950 -3.9458033 -6.6293585 -9.8115238
## [13] -3.2453047 -13.0205410 -4.7825523 -10.7126712 -13.3766017 -4.8749376
## [19] 1.0817535 1.2240558 -13.8441007 0.3954491 -4.8352020 -0.7766540
## [25] -3.0278188 -16.9286339 -3.7971107 -0.3820809 5.5137070 -5.7868252
## [31] 4.6410937 5.7645480 2.1694404 5.5117879 -14.0420511 6.2108247
## [37] 6.8589286 3.7200152 3.2310016 -10.3782165 -0.9952331 -11.8646028
## [43] 5.3300178 2.5262788 -39.3143246
##
## $dbz
## [1] -9.5445338 -9.9516071 -10.4835391 -10.8940401 -10.8991067 -10.3934633
## [7] -9.5467663 -8.6219866 -7.8072345 -7.1934660 -6.7976916 -6.5685554
## [13] -6.3774430 -6.0404364 -5.4233839 -4.5607251 -3.6260011 -2.7990140
## [19] -2.1908349 -1.8421559 -1.7384541 -1.8120718 -1.9298323 -1.8943331
## [25] -1.5187546 -0.7660644 0.2171424 1.2281916 2.1231629 2.8387644
## [31] 3.3707218 3.7473861 4.0052645 4.1670331 4.2285266 4.1620034
## [37] 3.9342417 3.5311982 2.9829486 2.3823143 1.8743874 1.5844664
## [43] 1.5205176 1.5702303 1.6044286We can see that accounting for a trend reducing the wandering behavior of the model.
aic5.wge(p.ns, p=0:10, q=0:2) # arima(10,1,0)## ---------WORKING... PLEASE WAIT...
##
##
## Five Smallest Values of aic## p q aic
## 31 10 0 -7.507572
## 33 10 2 -7.489151
## 32 10 1 -7.489061
## 29 9 1 -7.440917
## 7 2 0 -7.438540aic5.wge(p.ns, p=0:10, q=0:2, type='bic') #arima(2,1,0)## ---------WORKING... PLEASE WAIT...
##
##
## Five Smallest Values of bic## p q bic
## 7 2 0 -7.355213
## 2 0 1 -7.342572
## 5 1 1 -7.332269
## 3 0 2 -7.323260
## 8 2 1 -7.305147We can see that the AIC picks a 10,1,0 model while the bic picked the 2,1,0 model. We will move fowrad with the bic model as it is simpler to work with.
Next, we will apply the p.ns differenced data to estimate based on our ARIMA output.
estp.ns = est.arma.wge(p.ns,p=2,q=0)##
## Coefficients of Original polynomial:
## -0.7757 -0.3970
##
## Factor Roots Abs Recip System Freq
## 1+0.7757B+0.3970B^2 -0.9769+-1.2508i 0.6301 0.3555
##
## estp.ns## $phi
## [1] -0.7756653 -0.3969916
##
## $theta
## [1] 0
##
## $res
## [1] -1.188084e-03 -1.711727e-02 7.232104e-02 2.902631e-02 -6.655573e-02
## [6] 5.021585e-02 -1.639603e-02 2.655964e-02 4.919771e-03 2.732184e-03
## [11] 4.124697e-02 2.819740e-02 -2.541397e-02 -3.496573e-02 5.188099e-03
## [16] 3.014440e-02 1.906693e-02 4.619081e-02 -6.672304e-03 2.249080e-02
## [21] 3.146378e-02 2.796105e-02 1.269114e-02 -8.996794e-02 -1.199991e-03
## [26] -5.563065e-03 2.046206e-02 7.810951e-03 6.079474e-03 1.657065e-02
## [31] 1.840531e-02 -9.269985e-03 -2.589223e-02 -3.526198e-02 3.514342e-02
## [36] 3.205585e-03 1.062591e-02 -1.039418e-02 -1.344943e-02 -1.267713e-02
## [41] -8.889166e-03 -8.500580e-02 1.139347e-02 8.006710e-03 1.936007e-03
## [46] -1.440245e-02 -1.849550e-02 -3.420213e-03 -1.017148e-02 1.252846e-02
## [51] -4.874298e-03 2.923439e-03 -8.149795e-03 -1.396376e-02 -2.493752e-02
## [56] 6.033954e-03 5.753127e-04 2.999353e-03 -8.504751e-03 -7.274692e-03
## [61] -3.705642e-02 7.370685e-04 2.169728e-03 1.403475e-02 5.371820e-04
## [66] 9.400173e-04 -1.384365e-02 -3.873463e-03 -6.632044e-03 -9.360112e-04
## [71] 4.839194e-03 9.035012e-03 1.712204e-03 -2.895536e-03 -9.786221e-03
## [76] 4.788193e-03 1.523541e-02 2.022898e-03 -6.426013e-03 2.676855e-03
## [81] 5.512138e-03 -6.966347e-04 -3.872233e-03 -5.273892e-03 5.720479e-05
## [86] 2.900075e-03 2.697955e-03 -1.088815e-03 4.522679e-04 -2.359587e-04
##
## $avar
## [1] 0.000550212
##
## $aic
## [1] -7.43854
##
## $aicc
## [1] -6.411089
##
## $bic
## [1] -7.355213
##
## $se.phi
## [1] 0.009237752 0.009103055
##
## $se.theta
## [1] 0mean(p.ns)## [1] -0.001079281The preferred estimated model is (1 + 0.78B + 0.40B^2)(X - B) = a_t
fore.arma.wge(p.ns,phi=estp.ns$phi, n.ahead=10,plot=TRUE,limits=FALSE)
## $f
## [1] -0.0010944305 -0.0006846284 -0.0013793850 -0.0010031744 -0.0010191754
## [6] -0.0011561164 -0.0010435438 -0.0010764980 -0.0010956269 -0.0010677067
##
## $ll
## [1] -0.04745045 -0.05935122 -0.06080817 -0.06083296 -0.06154165 -0.06183377
## [7] -0.06172176 -0.06178521 -0.06182022 -0.06179311
##
## $ul
## [1] 0.04526159 0.05798196 0.05804940 0.05882662 0.05950330 0.05952154
## [7] 0.05963467 0.05963221 0.05962897 0.05965770
##
## $resid
## [1] 0.000000e+00 0.000000e+00 7.232104e-02 2.902631e-02 -6.655573e-02
## [6] 5.021585e-02 -1.639603e-02 2.655964e-02 4.919771e-03 2.732184e-03
## [11] 4.124697e-02 2.819740e-02 -2.541397e-02 -3.496573e-02 5.188099e-03
## [16] 3.014440e-02 1.906693e-02 4.619081e-02 -6.672304e-03 2.249080e-02
## [21] 3.146378e-02 2.796105e-02 1.269114e-02 -8.996794e-02 -1.199991e-03
## [26] -5.563065e-03 2.046206e-02 7.810951e-03 6.079474e-03 1.657065e-02
## [31] 1.840531e-02 -9.269985e-03 -2.589223e-02 -3.526198e-02 3.514342e-02
## [36] 3.205585e-03 1.062591e-02 -1.039418e-02 -1.344943e-02 -1.267713e-02
## [41] -8.889166e-03 -8.500580e-02 1.139347e-02 8.006710e-03 1.936007e-03
## [46] -1.440245e-02 -1.849550e-02 -3.420213e-03 -1.017148e-02 1.252846e-02
## [51] -4.874298e-03 2.923439e-03 -8.149795e-03 -1.396376e-02 -2.493752e-02
## [56] 6.033954e-03 5.753127e-04 2.999353e-03 -8.504751e-03 -7.274692e-03
## [61] -3.705642e-02 7.370685e-04 2.169728e-03 1.403475e-02 5.371820e-04
## [66] 9.400173e-04 -1.384365e-02 -3.873463e-03 -6.632044e-03 -9.360112e-04
## [71] 4.839194e-03 9.035012e-03 1.712204e-03 -2.895536e-03 -9.786221e-03
## [76] 4.788193e-03 1.523541e-02 2.022898e-03 -6.426013e-03 2.676855e-03
## [81] 5.512138e-03 -6.966347e-04 -3.872233e-03 -5.273892e-03 5.720479e-05
## [86] 2.900075e-03 2.697955e-03 -1.088815e-03 4.522679e-04 -2.359587e-04
##
## $wnv
## [1] 0.0005593712
##
## $se
## [1] 0.02365103 0.02993193 0.03032081 0.03052540 0.03087881 0.03095799
## [7] 0.03095827 0.03097383 0.03098194 0.03098235
##
## $psi
## [1] -0.775665268 0.204665039 0.149181010 -0.196964823 0.093555169
## [6] 0.005625879 -0.041504412 0.029960105 -0.006762111 -0.006648774forNS = fore.arma.wge(p.ns,phi=estp.ns$phi, n.ahead=10,lastn=TRUE,plot=TRUE,limits=FALSE)
Both forecasts accounting for the trend will oscillate back towards the mean, but since our model is differenced and it seems like trend is helping the forecast, the oscillating is minimal.
Now let’s pressure test the findings with an ASE score.
ASENS = mean((p.ns[(90-20+1):90] - forNS$f)^2)
ASENS## [1] 3.540763e-05ASE for data accounting for a trend ended up being much lower than a stationary model. However, when we run the same model on the data nonstationarized, we get a much higher ASE.
trainingSize = 70
horizon = 10
ASEHolder = numeric()
for( i in 1:(91-(trainingSize + horizon) + 1))
{
forecasts = fore.aruma.wge(p.ns[i:(i+(trainingSize-1))],phi = estp.ns$phi, theta = estp.ns$theta, s = 0, d = 1,n.ahead = horizon)
ASE = mean((p.ns[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
ASEHolder[i] = ASE
}
ASEHolder## [1] 5.211272e-05 4.878553e-05 6.328122e-05 4.747003e-05 5.176300e-05
## [6] 6.862660e-05 4.189509e-05 7.324802e-05 1.960326e-05 2.022338e-05
## [11] 7.861840e-06 NAhist(ASEHolder)
WindowedASE = mean(ASEHolder)
summary(ASEHolder)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 7.90e-06 3.11e-05 4.88e-05 4.50e-05 5.77e-05 7.32e-05 1WindowedASE## [1] NAOn the transformed realization, we see that we get a really low ASE, which isn’t totally fair since both are tracking back towards the mean and our realization is doing the same due to the nature of the covid trend.
Next we will look at a seasonal model. To do so, we will difference based on a week-long seasonality which was the most apparent in our EDA.
# Beginning the process of building a seasonal model (s=7)
y=artrans.wge(posper, phi.tr=c(0,0,0,0,0,0,1))
y_difTwice=artrans.wge(p.ns, phi.tr=c(0,0,0,0,0,0,1))
The seasonality forecast also reduces the wandering behavior of our model.
aic5.wge(y_difTwice, p = 0:12) # selects an aruma (10,1,0)s=7## ---------WORKING... PLEASE WAIT...
##
##
## Five Smallest Values of aic## p q aic
## 31 10 0 -7.208445
## 37 12 0 -7.195797
## 34 11 0 -7.193932
## 30 9 2 -7.188267
## 32 10 1 -7.188208aic5.wge(y_difTwice, p = 0:12, type='bic') #selects a aruma(9,1,0)s=7## ---------WORKING... PLEASE WAIT...
##
##
## Five Smallest Values of bic## p q bic
## 28 9 0 -6.888613
## 31 10 0 -6.887876
## 25 8 0 -6.881006
## 24 7 2 -6.862468
## 29 9 1 -6.858339The model we will be using to estimate our phis and thetas will be an ARUMA(9,1,0) with a seasonality of 7.
esty = est.arma.wge(y_difTwice,p=9, q=0)##
## Coefficients of Original polynomial:
## -0.9141 -0.6719 -0.3358 -0.2939 -0.2154 -0.1031 -0.6402 -0.5480 -0.2406
##
## Factor Roots Abs Recip System Freq
## 1-0.4277B+0.9114B^2 0.2346+-1.0208i 0.9547 0.2140
## 1+1.1155B+0.8946B^2 -0.6235+-0.8539i 0.9459 0.3504
## 1-1.5783B+0.8331B^2 0.9472+-0.5505i 0.9127 0.0838
## 1+0.9121B -1.0964 0.9121 0.5000
## 1+0.8924B+0.3882B^2 -1.1493+-1.1202i 0.6231 0.3770
##
## esty$phi## [1] -0.9141070 -0.6718829 -0.3357585 -0.2938873 -0.2154305 -0.1031342 -0.6401727
## [8] -0.5480305 -0.2405525esty$theta## [1] 0factor.wge(phi=esty$phi)##
## Coefficients of Original polynomial:
## -0.9141 -0.6719 -0.3358 -0.2939 -0.2154 -0.1031 -0.6402 -0.5480 -0.2406
##
## Factor Roots Abs Recip System Freq
## 1-0.4277B+0.9114B^2 0.2346+-1.0208i 0.9547 0.2140
## 1+1.1155B+0.8946B^2 -0.6235+-0.8539i 0.9459 0.3504
## 1-1.5783B+0.8331B^2 0.9472+-0.5505i 0.9127 0.0838
## 1+0.9121B -1.0964 0.9121 0.5000
## 1+0.8924B+0.3882B^2 -1.1493+-1.1202i 0.6231 0.3770
##
## The preferred estimated model is (1 + 0.91B + 0.67B^2 + 0.34B^3 + 0.29B^4 + 0.22B^5 + 0.10B^6 + 0.64B^7 + 0.55B^8 + 0.24B^9)(X - B)(X - B^7) = a_t
fore.aruma.wge(posper,phi=esty$phi, theta=esty$theta, d=1,s=7,n.ahead=20, lastn=FALSE, plot=TRUE)
## $f
## [1] 0.04752762 0.04112902 0.03497783 0.04110108 0.04001353 0.03598494
## [7] 0.03475024 0.03530077 0.03347623 0.03358888 0.03481040 0.03113134
## [13] 0.03089442 0.02907354 0.03272395 0.02856843 0.02382178 0.02823813
## [19] 0.02734893 0.02429649
##
## $ll
## [1] 0.0001412081 -0.0064318752 -0.0140316966 -0.0110851445 -0.0127460744
## [6] -0.0183224467 -0.0216299633 -0.0290883943 -0.0328487980 -0.0375576698
## [11] -0.0400061427 -0.0462122760 -0.0497711870 -0.0540502383 -0.0678083862
## [16] -0.0750312307 -0.0836352200 -0.0846741294 -0.0890216915 -0.0961214623
##
## $ul
## [1] 0.09491404 0.08868991 0.08398736 0.09328731 0.09277314 0.09029232
## [7] 0.09113044 0.09968993 0.09980127 0.10473544 0.10962693 0.10847496
## [13] 0.11156003 0.11219732 0.13325628 0.13216808 0.13127878 0.14115039
## [19] 0.14371954 0.14471445
##
## $resid
## [1] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
## [6] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
## [11] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
## [16] 0.0000000000 0.0000000000 -0.0290550639 0.0009558096 0.0324217115
## [21] 0.0066536103 0.0479074788 0.0296959084 0.0144988200 -0.0842028172
## [26] -0.0471520183 -0.0410674377 0.0004830700 0.0021912008 -0.0230942343
## [31] -0.0194938833 0.0344321550 -0.0342795751 -0.0038946757 -0.0448163019
## [36] 0.0107766274 -0.0188052962 -0.0103429763 0.0331524176 -0.0079832078
## [41] -0.0109909097 -0.0106329748 -0.1039270903 -0.0029619984 -0.0222427973
## [46] -0.0186943547 -0.0135162702 -0.0028947632 0.0229839332 0.0035453026
## [51] 0.0072344589 0.0065183604 0.0158567461 0.0120826919 0.0066569431
## [56] -0.0185424132 0.0605571443 -0.0059106098 -0.0056871036 -0.0109196432
## [61] 0.0044728606 -0.0221375576 0.0083349834 0.0034837235 0.0080729524
## [66] 0.0028317404 0.0063228737 -0.0028941642 0.0211503140 0.0122332862
## [71] 0.0020372420 0.0002439793 0.0087145528 0.0067026522 0.0069179202
## [76] 0.0152234055 0.0111595255 0.0146791506 -0.0032635838 -0.0062776721
## [81] 0.0026850683 0.0144444174 0.0061965052 -0.0016979753 -0.0087875952
## [86] -0.0050991172 -0.0046836139 -0.0022564577 -0.0021298770 0.0055776891
## [91] -0.0015605383
##
## $wnv
## [1] 0.0005845149
##
## $se
## [1] 0.02417674 0.02426576 0.02500486 0.02662563 0.02691817 0.02770785
## [7] 0.02876541 0.03285161 0.03383930 0.03629926 0.03817170 0.03946103
## [13] 0.04115592 0.04241009 0.05129201 0.05285697 0.05482500 0.05760830
## [19] 0.05937276 0.06143773
##
## $psi
## [1] 0.08589296 0.24960173 0.37836877 0.16370072 0.27166129 0.31965134
## [7] 0.65631523 0.33569577 0.54330619 0.48842712 0.41381491 0.48348871
## [13] 0.42344309 1.19324653 0.52805793 0.60212937 0.73169383 0.59424681
## [19] 0.65329775 0.70850934
##
## $ptot
## [1] 17
##
## $phitot
## [1] 0.085892960 0.242224128 0.336124388 0.041871225 0.078456830
## [6] 0.112296292 0.462961524 0.006249217 0.065253894 -0.095571933
## [11] -0.041871225 -0.078456830 -0.112296292 0.537038476 -0.092142177
## [16] -0.307478022 -0.240552455forS7 = fore.aruma.wge(posper,phi=esty$phi, theta=esty$theta, d=1,s=7,n.ahead=20, lastn=TRUE, plot=TRUE)
Now, let’s check our ASE for the Seasonal Trend model.
ASEs7 = mean((posper[(90-20+1):90] - forS7$f)^2)
ASEs7## [1] 0.0007639599Here, we did pretty well with our ASE score on a 20 day lookback.
trainingSize = 70
horizon = 10
ASEHolder = numeric()
for( i in 1:(81-(trainingSize + horizon) + 1))
{
forecasts = fore.aruma.wge(posper[i:(i+(trainingSize-1))],phi = esty$phi, theta = esty$theta, s = 7, d = 1,n.ahead = horizon)
ASE = mean((posper[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
ASEHolder[i] = ASE
}
ASEHolder## [1] 0.0002189011 0.0002497114hist(ASEHolder)
WindowedASE = mean(ASEHolder)
summary(ASEHolder)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0002189 0.0002266 0.0002343 0.0002343 0.0002420 0.0002497WindowedASE## [1] 0.0002343063So in Summary, our ASE scores are as follows.
print('So, In Summary, the ASE Values Are:')## [1] "So, In Summary, the ASE Values Are:"print(paste('Arma(2,2):', ASE22))## [1] "Arma(2,2): 0.000259881145874531"print(paste('Arima(2,1,0):', ASENS))## [1] "Arima(2,1,0): 3.54076321872339e-05"print(paste('Aruma(9,1,0)s=7:', ASEs7))## [1] "Aruma(9,1,0)s=7: 0.000763959941806105"The best performing model with the US data is the ARIMA(2,1,0). Interestingly, while we see what visually appears to be a seasonal trend, accounting for it does not produce the best mode.
Var Forecast
For our multivariate forecasts, we will use the VAR model to predict positive percentage in the US and AZ data.
To get started, let’s look into seeing if there’s a lag between other variables and the ccf.
ccf(posper, us2$death[36:182]) # lag is 1
ccf(posper, us2$hospitalizedCurrently[56:182]) # lag is 10
ccf(posper, us2$inIcuCurrently[65:182]) #lag is 3
ccf(posper, us2$onVentilatorCurrently[64:182]) # lag is -12
ccf(posper, us2$deathIncrease) #lag is -11
ccf(posper, us2$totalTestResultsIncrease) #lag is -3
Var Forecast of Stationary data.
Looking at the CCFs of positive percentage with other key variables, we can see that almost all of them have some sort of lagged correlation with the residual variables. So when we run a multivariate forecast, we can account for these lags to strengthen the forecast.
us2[is.na(us2)] <- 0library(vars)## Loading required package: MASS## Loading required package: strucchange## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: sandwich## Loading required package: urca## Loading required package: lmtestus2$posper = posper
#us2$posper[1:172], us2$death[1:172], us2$hospitalizedCurrently[1:172], us2$inIcuCurrently[1:172], us2$onVentilatorCurrently[1:172], us2$deathIncrease[1:172], us2$totalTestResultsIncrease[1:172], us2$negativeIncrease[1:172], us2$positiveIncrease[1:172]
BSVar1 = VARselect(cbind(us2$posper[92:172], us2$totalTestResultsIncrease[92:172], us2$positiveIncrease[92:172]), type = 'both', lag.max = 5)
BSVar1 # AIC was 3.02e01 for p =5## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 5 3 2 5
##
## $criteria
## 1 2 3 4 5
## AIC(n) 3.074682e+01 3.040404e+01 3.027063e+01 3.034676e+01 3.024378e+01
## HQ(n) 3.093067e+01 3.069819e+01 3.067509e+01 3.086152e+01 3.086885e+01
## SC(n) 3.120684e+01 3.114006e+01 3.128266e+01 3.163480e+01 3.180782e+01
## FPE(n) 2.256440e+13 1.604456e+13 1.409404e+13 1.530995e+13 1.395501e+13Running a var that accounts for the lags between positive percentage, deaths, numbers of people in the ICU, the number of hospitalized patients, number of people on ventilators, number of positives, number of negatives and total test results, we can see that accounting for all the lags the preferred AIC is a p of 10 at an AIC of 3.58e+01.
Now let’s look to see if we can predict the forecast on the var.
#ustest <- data.frame(posper, us2$death, us2$hospitalizedCurrently, us2$inIcuCurrently, us2$onVentilatorCurrently, us2$totalTestResultsIncrease, us2$negativeIncrease, us2$positiveIncrease)
#ustest = cbind(us2$posper[40:140], us2$death[40:140], us2$hospitalizedCurrently[40:140], us2$inIcuCurrently[40:140], us2$onVentilatorCurrently[40:140], us2$totalTestResultsIncrease[40:140], us2$negativeIncrease[40:140], us2$positiveIncrease[40:140])
PosperVAR = VAR(cbind(us2$posper[92:172], us2$totalTestResultsIncrease[92:172], us2$positiveIncrease[92:172]), type = "const", p = 5)## Warning in VAR(cbind(us2$posper[92:172], us2$totalTestResultsIncrease[92:172], : No column names supplied in y, using: y1, y2, y3 , instead.preds10stat=predict(PosperVAR,n.ahead=10)
preds20stat=predict(PosperVAR,n.ahead=20)In order to be able to plot our findings, we will need to prepare our data by appending our forecasts on our original data.
startingPoints10 = us2$posper[163:172]
posperForecasts10 = preds10stat$fcst$y1[,1:4] + startingPoints10
startingPoints20 = us2$posper[153:172]
posperForecasts20 = preds20stat$fcst$y1[,1:4] + startingPoints20Now let’s plot our findings using our forecasts overlayed on top of our original data.
#dev.off()
library(RColorBrewer)
#fanchart(preds, colors = brewer.pal(n = 8, name= 'Blues'))
plot(seq(92,182,1), us2$posper[92:182], type = "l",xlim = c(92,185), ylim = c(0,0.6), ylab = "Positive Percentage", main = "10 Day Positive Percentage Forecast")
lines(seq(173,182,1), as.data.frame(posperForecasts10)$fcst, type = "l", col = "red")
plot(seq(92,182,1), us2$posper[92:182], type = "l",xlim = c(92,185), ylim = c(0,0.6), ylab = "Positive Percentage", main = "10 Day Positive Percentage Forecast")
lines(seq(163,182,1), as.data.frame(posperForecasts20)$fcst, type = "l", col = "red")
Using the Var forecast including death, hospitalized patients, patients on ventilators and patients in ICU, we can see that the prediction over emphasize the spike that would happen between observations 172 - 182. Lets check the ASE to see how we did.
ASEvarstat = mean((us2$posper[173:182] - posperForecasts10[,1])^2)
ASEvarstat #0.0804## [1] 0.002822875ASEvarstat20 = mean((us2$posper[163:182] - posperForecasts20[,1])^2)
ASEvarstat20 #0.0804## [1] 0.003291832Now lets do a var assuming a trend model (which ended up being the best in a univariate study)
Var forecast with trend
pospertrend = artrans.wge(us2$posper, 1)
testresultstrends = artrans.wge(us2$totalTestResultsIncrease, 1)
positivetesttrend = artrans.wge(us2$positiveIncrease, 1)
Now that we have differenced variables, lets see which gives us the best AIC score
BSVar1 = VARselect(cbind(pospertrend[92:172], testresultstrends[92:172], positivetesttrend[92:172]), type = 'both', lag.max = 10)
BSVar1 # AIC was 3.03e01 for p =5## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 5 2 1 5
##
## $criteria
## 1 2 3 4 5
## AIC(n) 3.066764e+01 3.052586e+01 3.046363e+01 3.037874e+01 3.024397e+01
## HQ(n) 3.085774e+01 3.083002e+01 3.088185e+01 3.091102e+01 3.089030e+01
## SC(n) 3.114567e+01 3.129071e+01 3.151530e+01 3.171723e+01 3.186927e+01
## FPE(n) 2.084931e+13 1.813284e+13 1.711848e+13 1.585404e+13 1.403265e+13
## 6 7 8 9 10
## AIC(n) 3.034189e+01 3.027760e+01 3.024597e+01 3.035901e+01 3.036134e+01
## HQ(n) 3.110228e+01 3.115205e+01 3.123448e+01 3.146158e+01 3.157797e+01
## SC(n) 3.225401e+01 3.247655e+01 3.273174e+01 3.313159e+01 3.342074e+01
## FPE(n) 1.576556e+13 1.516820e+13 1.520943e+13 1.781053e+13 1.890829e+13Like the earlier model, we end up with an AIC preference of 5 with the AIC of 3.02.
PosperVAR = VAR(cbind(pospertrend[92:172], testresultstrends[92:172], positivetesttrend[92:172]), type = "const", p = 5)## Warning in VAR(cbind(pospertrend[92:172], testresultstrends[92:172], positivetesttrend[92:172]), : No column names supplied in y, using: y1, y2, y3 , instead.preds10=predict(PosperVAR,n.ahead=10)
preds20=predict(PosperVAR,n.ahead=20)startingPoints10 = pospertrend[130:139]
posperForecasts10 = preds10$fcst$y1[,1:3] + startingPoints10
startingPoints20 = pospertrend[120:139]
posperForecasts20 = preds20$fcst$y1[,1:3] + startingPoints20Then we will prepare our data so that we can overlay our forecast on top of the actual data.
#dev.off()
library(RColorBrewer)
#fanchart(preds, colors = brewer.pal(n = 8, name= 'Blues'))
plot(seq(1,139,1), pospertrend[1:139], type = "l",xlim = c(1,145), ylim = c(-0.5,0.5), ylab = "Positive Percentage", main = "10 Day Positive Percentage Forecast")
lines(seq(130,139,1), as.data.frame(posperForecasts10)$fcst, type = "l", col = "red")
plot(seq(1,139,1), pospertrend[1:139], type = "l",xlim = c(1,145), ylim = c(-0.5,0.5), ylab = "Positive Percentage", main = "10 Day Positive Percentage Forecast")
lines(seq(120,139,1), as.data.frame(posperForecasts20)$fcst, type = "l", col = "red")
Here in our model, we can see that our forecast over emphasized the spikes that appeared earlier in our data when, in turn, didn’t appear in the actual data.
ASEvartrend = mean((pospertrend[130:139] - posperForecasts10[,1])^2)
ASEvartrend #0.0036## [1] 0.0004515429ASEvartrend20 = mean((pospertrend[120:139] - posperForecasts20[,1])^2)
ASEvartrend20 #0.0036## [1] 0.0002274616We can see with a trend model, we get a better ASE than than the stationary model we ran previously.
Now lets try a seasonal model accounting for 7 days on top of the trend data we previously ran.
pospertrend2 = artrans.wge(pospertrend, c(rep(0,6),1))
testresultstrends2 = artrans.wge(testresultstrends, c(rep(0,6),1))
positivetesttrend2 = artrans.wge(positivetesttrend, c(rep(0,6),1))
For our airplane model, we can see that we are differencing our trend data based on a seasonality factor of 7 days, which we saw in the previous EDA. Here, we can see that our low pass filter created more emphasis on the higher frequency periods later in the dataframe.
BSVar1 = VARselect(cbind(pospertrend2[92:172], testresultstrends2[92:172], positivetesttrend2[92:172]), type = 'both', lag.max = 10)
BSVar1 # AIC was 3.09e01 for p =9## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 9 2 1 8
##
## $criteria
## 1 2 3 4 5
## AIC(n) 3.153286e+01 3.130455e+01 3.133644e+01 3.132847e+01 3.142546e+01
## HQ(n) 3.172296e+01 3.160871e+01 3.175466e+01 3.186075e+01 3.207179e+01
## SC(n) 3.201089e+01 3.206940e+01 3.238811e+01 3.266696e+01 3.305076e+01
## FPE(n) 4.952825e+13 3.950455e+13 4.097527e+13 4.098284e+13 4.573559e+13
## 6 7 8 9 10
## AIC(n) 3.117114e+01 3.106300e+01 3.101419e+01 3.099480e+01 3.117415e+01
## HQ(n) 3.193153e+01 3.193745e+01 3.200270e+01 3.209736e+01 3.239078e+01
## SC(n) 3.308327e+01 3.326195e+01 3.349995e+01 3.376738e+01 3.423355e+01
## FPE(n) 3.612849e+13 3.326823e+13 3.279032e+13 3.363519e+13 4.262360e+13PosperVAR = VAR(cbind(pospertrend2[92:172], testresultstrends[92:172], positivetesttrend[92:172]), type = "const", p = 9)## Warning in VAR(cbind(pospertrend2[92:172], testresultstrends[92:172], positivetesttrend[92:172]), : No column names supplied in y, using: y1, y2, y3 , instead.preds=predict(PosperVAR,n.ahead=10)After differencing, we can see that our Varselect is selecting an AIC of 9 at the level of 3.09e+01. And to do so, we made sure our positive percentage only applies to the data that related best to them.
startingPoints10 = pospertrend2[163:172]
posperForecasts10 = preds10$fcst$y1[,1:3] + startingPoints10
startingPoints20 = pospertrend2[153:172]
posperForecasts20 = preds20$fcst$y1[,1:3] + startingPoints20Now that we ran our forecast, we will need to prepare our data so we can visualize performance.
#dev.off()
library(RColorBrewer)
#fanchart(preds, colors = brewer.pal(n = 8, name= 'Blues'))
plot(seq(92,174,1), pospertrend2[92:174], type = "l",xlim = c(92,180), ylim = c(-0.5,0.5), ylab = "Positive Percentage", main = "10 Day Positive Percentage Forecast")
lines(seq(165,174,1), as.data.frame(posperForecasts10)$fcst, type = "l", col = "red")
plot(seq(92,174,1), pospertrend2[92:174], type = "l",xlim = c(92,180), ylim = c(-0.5,0.5), ylab = "Positive Percentage", main = "10 Day Positive Percentage Forecast")
lines(seq(155,174,1), as.data.frame(posperForecasts20)$fcst, type = "l", col = "red")
Now that we have run our Var, we can see that our 10 and 20 day lookback window is looking strong, however with the extra difference, we can see it’s over emphasizing the seasonality spikes.
ASEvarairplane = mean((pospertrend2[123:132] - posperForecasts10[,1])^2)
ASEvarairplane #0.0021## [1] 0.004647224ASEvarairplane20 = mean((pospertrend2[113:132] - posperForecasts20[,1])^2)
ASEvarairplane20 #0.0021## [1] 0.003393451We can see that our ASE score for the short and long term window are 0.014 and 0.010 respectively for 10 and 20 days.
print('So, In Summary, the ASE Values Are:')## [1] "So, In Summary, the ASE Values Are:"print(paste('Arma(2,2):', ASE22))## [1] "Arma(2,2): 0.000259881145874531"print(paste('Arima(2,1,0):', ASENS))## [1] "Arima(2,1,0): 3.54076321872339e-05"print(paste('Aruma(9,1,0)s=7:', ASEs7))## [1] "Aruma(9,1,0)s=7: 0.000763959941806105"print(paste('Var Stationary:', ASEvarstat))## [1] "Var Stationary: 0.00282287523520302"print(paste('Var Trend:', ASEvartrend))## [1] "Var Trend: 0.000451542857873454"print(paste('Var Airplane:', ASEvarairplane))## [1] "Var Airplane: 0.00464722379532156"Here we can see the best performing models are continuing to be focusing on mitigating a wandering trend rather than a stationary and seasonal/arplane models.
Next we will use mlp and neural nets to see if we can improve upon performance.
MLP
We will break the data set into 10 and 20 day lookback windows and difference our data based on trend into an ARIMA model.
We will then do the same for positive percentage and then test results and daily positives.
library(nnfor)## Loading required package: forecast## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(forecast)
library(vars)
library(tswge)
us2Small = us2[92:172,]
us2Large = us2[92:162,]
pospersmall = artrans.wge(us2Small$posper, phi.tr=(1))
testresultssmall = artrans.wge(us2Small$totalTestResultsIncrease, phi.tr=(1))
positivesmall = artrans.wge(us2Small$positiveIncrease, phi.tr=(1))
posperlarge = artrans.wge(us2Large$posper, phi.tr=(1))
testresultslarge = artrans.wge(us2Large$totalTestResultsIncrease, phi.tr=(1))
positivelarge = artrans.wge(us2Large$positiveIncrease, phi.tr=(1))
us2SmallDF = data.frame(test = ts(testresultssmall, frequency = 1), positive = ts(positivesmall, frequency = 1))
us2LargeDF = data.frame(test = ts(testresultslarge, frequency = 1), positive = ts(positivelarge, frequency = 1))From here, we will fit it into a neural net model.
fit.mlp.small = mlp(ts(pospersmall, frequency = 1), allow.det.season = FALSE, reps = 20, comb = 'mean', xreg = us2SmallDF)
#fit.mlp.small = mlp(ts(us2Small$deathIncrease, frequency = 1), allow.det.season = FALSE, reps = 20, comb = 'mean', xreg = us2SmallDF)
fit.mlp.small## MLP fit with 5 hidden nodes and 20 repetitions.
## Univariate lags: (1,2)
## 2 regressors included.
## - Regressor 1 lags: (3,4)
## - Regressor 2 lags: (1,2,4)
## Forecast combined using the mean operator.
## MSE: 0.plot(fit.mlp.small)
fit.mlp.large = mlp(ts(posperlarge, frequency = 1), allow.det.season = FALSE, reps = 20, comb = 'mean', xreg = us2LargeDF)
#fit.mlp.large = mlp(ts(us2Large$deathIncrease, frequency = 1), allow.det.season = FALSE, reps = 20, comb = 'mean', xreg = us2LargeDF)
fit.mlp.large## MLP fit with 5 hidden nodes and 20 repetitions.
## Univariate lags: (1,2)
## 2 regressors included.
## - Regressor 1 lags: (3,4)
## - Regressor 2 lags: (1,4)
## Forecast combined using the mean operator.
## MSE: 0.plot(fit.mlp.large)
The MLP output gave us 5 notes and 20 repititions, with two regressors included.
From here we will create a dataframe using the key features of the regressors, and then forecast against it using a horizon of 10 days.
testdif = artrans.wge(us2$totalTestResults, phi.tr=(1))
posdif = artrans.wge(us2$positiveIncrease, phi.tr=(1))
USDF = data.frame(test = ts(testdif), positive = ts(posdif))
fore.mlp10 = forecast(fit.mlp.small, h = 10, xreg = USDF)
fore.mlp20 = forecast(fit.mlp.large, h = 20, xreg = USDF)
plot(fore.mlp10)
plot(fore.mlp20)
We can see with our forecast model our expected values for a 10 and 20 day window. Our MLP continues to emphasize the spikes which accounts for the most variability with our forecast.
Now let’s check the ASE of our short term and long term forecast.
ASEmlp10 = mean((pospersmall[71:80] - fore.mlp10$mean)^2)
ASEmlp10## [1] 8.695629e-05ASEmlp20 = mean((posperlarge[51:70] - fore.mlp20$mean)^2)
ASEmlp20## [1] 0.0003007022Our 10 and 20 day forecast windows ASE’s for the MLP were 0.00018 and 0.00051 respectively.
Next we will combine our stationary forecast model with our MLP forecast model to see if we can improve our performance.
Ensemble Model
Here our ensemble will combine the forecast and stationary Univariate model. In doing so, we will look at a 10 day and 20 day lookback window.
ensemble10 = (preds10stat$fcst$y1[,1] + fore.mlp10$mean)/2
ensemble20 = (preds20stat$fcst$y1[,1] + fore.mlp20$mean)/2Now that we performed our ensemble, we will combine with our actual data frame in a 10 and 20 day lookback window.
plot(seq(92,182,1), us2$posper[92:182], type = "l",xlim = c(92,190), ylim = c(0,0.3), ylab = "Positive Percentage", main = "10 Day Positive Percentage Forecast Ensemble")
lines(seq(173,182,1), ensemble10, type = "l", col = "red")
plot(seq(92,182,1), us2$posper[92:182], type = "l",xlim = c(92,190), ylim = c(0,0.3), ylab = "Positive Percentage", main = "20 Day Positive Percentage Forecast Ensemble")
lines(seq(163,182,1), ensemble20, type = "l", col = "red")
Both our MLP and Ensemble forecast appears to be predicting our actuals to be closer to 0 than our other models.
Ensemble ASE
ASE10 = mean((us2$posper[173:182] - ensemble10)^2)
ASE10## [1] 0.0007792512ASE20 = mean((us2$posper[163:182] - ensemble20)^2)
ASE20## [1] 0.00110013With our respective ASE,s we can see both are hoving around 0.00089 in our 10 and 20 day lookback window.
print('So, In Summary, the ASE For our US Analysis Are:')## [1] "So, In Summary, the ASE For our US Analysis Are:"print(paste('Arma(2,2):', ASE22))## [1] "Arma(2,2): 0.000259881145874531"print(paste('Arima(2,1,0):', ASENS))## [1] "Arima(2,1,0): 3.54076321872339e-05"print(paste('Aruma(9,1,0)s=7:', ASEs7))## [1] "Aruma(9,1,0)s=7: 0.000763959941806105"print(paste('Var Stationary:', ASEvarstat))## [1] "Var Stationary: 0.00282287523520302"print(paste('Var Trend:', ASEvartrend))## [1] "Var Trend: 0.000451542857873454"print(paste('Var Airplane:', ASEvarairplane))## [1] "Var Airplane: 0.00464722379532156"print(paste('MLP:', ASEmlp10))## [1] "MLP: 8.69562904868902e-05"print(paste('Ensemble:', ASE10))## [1] "Ensemble: 0.000779251232977854"Looking across all our single ASE scores, we can see that our Arima (2,1,0) model had the best performing ASE score for the univariate data. When introducing multivariate data, we can see that our MLP model had the best performing ASE score.
Overall, it’s reassuring to see that when we forecast US data, our estimations are relatively close to how they actually performed.
Now Let’s look at the AZ data.
Arizona Data
plotts.sample.wge(as2$date)
## $autplt
## [1] 1.0000000 0.9785482 0.9569853 0.9353119 0.9135282 0.8916348 0.8696319
## [8] 0.8475201 0.8252996 0.8029710 0.7805345 0.7579907 0.7353399 0.7125825
## [15] 0.6897189 0.6667495 0.6436747 0.6204949 0.5972105 0.5738219 0.5503296
## [22] 0.5267338 0.5091127 0.4914171 0.4736474 0.4558039
##
## $freq
## [1] 0.007142857 0.014285714 0.021428571 0.028571429 0.035714286 0.042857143
## [7] 0.050000000 0.057142857 0.064285714 0.071428571 0.078571429 0.085714286
## [13] 0.092857143 0.100000000 0.107142857 0.114285714 0.121428571 0.128571429
## [19] 0.135714286 0.142857143 0.150000000 0.157142857 0.164285714 0.171428571
## [25] 0.178571429 0.185714286 0.192857143 0.200000000 0.207142857 0.214285714
## [31] 0.221428571 0.228571429 0.235714286 0.242857143 0.250000000 0.257142857
## [37] 0.264285714 0.271428571 0.278571429 0.285714286 0.292857143 0.300000000
## [43] 0.307142857 0.314285714 0.321428571 0.328571429 0.335714286 0.342857143
## [49] 0.350000000 0.357142857 0.364285714 0.371428571 0.378571429 0.385714286
## [55] 0.392857143 0.400000000 0.407142857 0.414285714 0.421428571 0.428571429
## [61] 0.435714286 0.442857143 0.450000000 0.457142857 0.464285714 0.471428571
## [67] 0.478571429 0.485714286 0.492857143 0.500000000
##
## $db
## [1] 16.263652 10.447498 7.087719 6.037123 -1.338094 -0.821000
## [7] -1.491848 -2.551596 -1.589645 -4.224602 -4.792586 -5.754813
## [13] -5.617517 -9.195804 -8.305557 -7.460668 -8.365526 -5.875778
## [19] -14.438691 -12.914235 -10.715221 -12.041493 -7.331212 -12.340935
## [25] -13.668936 -12.786866 -13.675082 -10.852331 -12.757535 -13.337064
## [31] -13.711836 -11.607137 -14.915976 -15.859237 -14.537489 -16.392025
## [37] -11.068526 -14.258853 -18.178917 -15.721681 -20.987353 -12.948443
## [43] -13.444602 -16.884671 -15.532006 -18.105579 -16.180284 -15.002900
## [49] -16.547571 -16.205770 -15.286917 -15.843107 -17.792538 -16.931200
## [55] -19.739346 -15.750383 -14.155637 -18.099023 -15.918000 -24.090310
## [61] -18.864051 -14.511381 -17.458014 -14.978069 -19.518007 -19.202924
## [67] -16.526817 -17.719266 -16.029290 -17.568569
##
## $dbz
## [1] 11.6577443 11.1469541 10.2886322 9.0734157 7.4917493 5.5414489
## [7] 3.2462956 0.6965271 -1.8889612 -4.1597523 -5.8526924 -7.0302261
## [13] -7.8946487 -8.5214835 -8.9080794 -9.1353881 -9.3733371 -9.7681011
## [19] -10.3707515 -11.1309757 -11.9263878 -12.6276042 -13.1763650 -13.6003020
## [25] -13.9484078 -14.2346451 -14.4478298 -14.5979355 -14.7358785 -14.9252887
## [31] -15.2009089 -15.5491698 -15.9190876 -16.2553695 -16.5316940 -16.7564471
## [37] -16.9491245 -17.1166409 -17.2531495 -17.3582067 -17.4499562 -17.5579248
## [43] -17.7030528 -17.8834531 -18.0765027 -18.2548639 -18.4039045 -18.5264273
## [49] -18.6326077 -18.7274035 -18.8079515 -18.8715046 -18.9238895 -18.9788173
## [55] -19.0484269 -19.1340115 -19.2250058 -19.3069340 -19.3716478 -19.4212422
## [61] -19.4630337 -19.5014494 -19.5347203 -19.5586937 -19.5730905 -19.5834622
## [67] -19.5967003 -19.6145041 -19.6311708 -19.6380863Here we can see that the date is in order which works for us in our realization.
plotts.sample.wge(as2$death[18:140])
## $autplt
## [1] 1.0000000 0.9647762 0.9318906 0.8983310 0.8641482 0.8328016 0.8028270
## [8] 0.7731511 0.7447694 0.7178591 0.6901319 0.6638545 0.6385270 0.6135843
## [15] 0.5898725 0.5661129 0.5448956 0.5229471 0.5003239 0.4773914 0.4546995
## [22] 0.4322101 0.4113310 0.3909803 0.3700409 0.3486550
##
## $freq
## [1] 0.008130081 0.016260163 0.024390244 0.032520325 0.040650407 0.048780488
## [7] 0.056910569 0.065040650 0.073170732 0.081300813 0.089430894 0.097560976
## [13] 0.105691057 0.113821138 0.121951220 0.130081301 0.138211382 0.146341463
## [19] 0.154471545 0.162601626 0.170731707 0.178861789 0.186991870 0.195121951
## [25] 0.203252033 0.211382114 0.219512195 0.227642276 0.235772358 0.243902439
## [31] 0.252032520 0.260162602 0.268292683 0.276422764 0.284552846 0.292682927
## [37] 0.300813008 0.308943089 0.317073171 0.325203252 0.333333333 0.341463415
## [43] 0.349593496 0.357723577 0.365853659 0.373983740 0.382113821 0.390243902
## [49] 0.398373984 0.406504065 0.414634146 0.422764228 0.430894309 0.439024390
## [55] 0.447154472 0.455284553 0.463414634 0.471544715 0.479674797 0.487804878
## [61] 0.495934959
##
## $db
## [1] 15.2651906 10.0014146 6.7008124 4.6478363 2.3725061 1.5784605
## [7] 0.1371372 -1.3749032 -2.3463539 -3.3103015 -4.0405719 -4.6738272
## [13] -5.3087895 -5.6679032 -6.2283494 -7.2127906 -7.3247616 -10.0919656
## [19] -8.7311147 -8.6777284 -9.0638777 -9.8972503 -10.8831717 -10.7020546
## [25] -11.3219069 -11.3588454 -11.5240570 -11.3390221 -12.1017645 -12.3781648
## [31] -12.6713416 -12.5390052 -13.3507635 -13.6140396 -14.1059594 -12.8640286
## [37] -12.9323299 -13.2556342 -13.3798885 -13.6103529 -13.8192301 -13.2266590
## [43] -13.6072158 -13.6944964 -13.7015604 -14.0941831 -14.5170015 -14.9947461
## [49] -14.9921989 -14.4377870 -15.0100001 -14.3210643 -14.4181132 -14.3888333
## [55] -14.5948606 -14.7866629 -15.5427642 -15.3883655 -15.4368991 -14.9470916
## [61] -14.8736730
##
## $dbz
## [1] 11.205587 10.627648 9.660324 8.301395 6.558657 4.469307
## [7] 2.137120 -0.222676 -2.308900 -3.924407 -5.125284 -6.060033
## [13] -6.785767 -7.316171 -7.735285 -8.183736 -8.762627 -9.477068
## [19] -10.244741 -10.954544 -11.544113 -12.025933 -12.437233 -12.787848
## [25] -13.067539 -13.288117 -13.497953 -13.750598 -14.065865 -14.418751
## [31] -14.759601 -15.046795 -15.264451 -15.414158 -15.500128 -15.529483
## [37] -15.522424 -15.512198 -15.529992 -15.589361 -15.684091 -15.799733
## [43] -15.927630 -16.068810 -16.225380 -16.389512 -16.541923 -16.662378
## [49] -16.743086 -16.792809 -16.829254 -16.868640 -16.920874 -16.990127
## [55] -17.075668 -17.170383 -17.260029 -17.327939 -17.364696 -17.374964
## [61] -17.373698For the case of deaths in Arizona, we can see that the model is increasing with the passing day which indicates that the death variable is a cumulative figure.
plotts.sample.wge(as2$positive)
## $autplt
## [1] 1.0000000 0.9681157 0.9364178 0.9039511 0.8711208 0.8381524 0.8056101
## [8] 0.7732416 0.7410177 0.7095172 0.6773215 0.6450589 0.6129986 0.5817596
## [15] 0.5513033 0.5213612 0.4921062 0.4634159 0.4354319 0.4077758 0.3813582
## [22] 0.3556582 0.3314692 0.3088227 0.2857131 0.2637792
##
## $freq
## [1] 0.007142857 0.014285714 0.021428571 0.028571429 0.035714286 0.042857143
## [7] 0.050000000 0.057142857 0.064285714 0.071428571 0.078571429 0.085714286
## [13] 0.092857143 0.100000000 0.107142857 0.114285714 0.121428571 0.128571429
## [19] 0.135714286 0.142857143 0.150000000 0.157142857 0.164285714 0.171428571
## [25] 0.178571429 0.185714286 0.192857143 0.200000000 0.207142857 0.214285714
## [31] 0.221428571 0.228571429 0.235714286 0.242857143 0.250000000 0.257142857
## [37] 0.264285714 0.271428571 0.278571429 0.285714286 0.292857143 0.300000000
## [43] 0.307142857 0.314285714 0.321428571 0.328571429 0.335714286 0.342857143
## [49] 0.350000000 0.357142857 0.364285714 0.371428571 0.378571429 0.385714286
## [55] 0.392857143 0.400000000 0.407142857 0.414285714 0.421428571 0.428571429
## [61] 0.435714286 0.442857143 0.450000000 0.457142857 0.464285714 0.471428571
## [67] 0.478571429 0.485714286 0.492857143 0.500000000
##
## $db
## [1] 15.40072875 11.84123324 8.14561086 4.97379521 2.80360046
## [6] 1.14335542 -0.02900222 -1.42461280 -2.43820173 -3.04581229
## [11] -3.80792669 -4.77316921 -5.52817564 -5.97223495 -6.40516131
## [16] -6.99158367 -7.59371107 -8.14153208 -8.50504366 -9.31828911
## [21] -10.09435302 -10.20877443 -10.22479227 -10.52165808 -10.81330872
## [26] -11.21516951 -11.62370242 -11.85128549 -11.73267761 -11.91283903
## [31] -12.33064733 -12.68759392 -13.16192728 -13.56178986 -13.37885323
## [36] -13.38559714 -13.64455375 -13.97889567 -14.35480611 -14.23683318
## [41] -13.92973555 -14.20735368 -14.51782639 -14.59027297 -14.83253177
## [46] -14.89163927 -15.12555130 -14.98019114 -14.87387781 -15.29222516
## [51] -15.73334474 -15.82731091 -15.96373193 -15.77413118 -15.40370988
## [56] -15.70425798 -16.11354839 -16.31835800 -16.25404623 -15.70251541
## [61] -15.47017803 -15.85315043 -16.09666670 -16.03925082 -16.10038959
## [66] -16.20032127 -16.20249575 -16.18266196 -16.23887770 -16.31161943
##
## $dbz
## [1] 11.3885341 10.9167836 10.1275618 9.0181834 7.5895603 5.8540199
## [7] 3.8507881 1.6718887 -0.5118262 -2.4717288 -4.0389303 -5.2139903
## [13] -6.1012862 -6.7801140 -7.2912677 -7.6951221 -8.0868993 -8.5547980
## [19] -9.1371627 -9.8095205 -10.5003600 -11.1296272 -11.6514679 -12.0673390
## [25] -12.4012188 -12.6714390 -12.8885170 -13.0706363 -13.2507340 -13.4644627
## [31] -13.7309173 -14.0423078 -14.3690687 -14.6767580 -14.9434209 -15.1649742
## [37] -15.3471071 -15.4951387 -15.6125213 -15.7067395 -15.7930728 -15.8901065
## [43] -16.0103599 -16.1540295 -16.3106027 -16.4667640 -16.6143758 -16.7525548
## [49] -16.8833541 -17.0062958 -17.1172316 -17.2120880 -17.2911312 -17.3590809
## [55] -17.4208402 -17.4770615 -17.5238339 -17.5568610 -17.5763954 -17.5885319
## [61] -17.6018757 -17.6226638 -17.6522229 -17.6880840 -17.7269441 -17.7665343
## [67] -17.8049863 -17.8389578 -17.8630652 -17.8718748plotts.sample.wge(as2$negative)
## $autplt
## [1] 1.0000000 0.9762416 0.9525335 0.9284555 0.9045013 0.8800691 0.8557717
## [8] 0.8314940 0.8073024 0.7836117 0.7595330 0.7354845 0.7116705 0.6881240
## [15] 0.6644983 0.6408512 0.6170985 0.5934766 0.5701367 0.5468273 0.5239327
## [22] 0.5010465 0.4785424 0.4567996 0.4347440 0.4130348
##
## $freq
## [1] 0.007142857 0.014285714 0.021428571 0.028571429 0.035714286 0.042857143
## [7] 0.050000000 0.057142857 0.064285714 0.071428571 0.078571429 0.085714286
## [13] 0.092857143 0.100000000 0.107142857 0.114285714 0.121428571 0.128571429
## [19] 0.135714286 0.142857143 0.150000000 0.157142857 0.164285714 0.171428571
## [25] 0.178571429 0.185714286 0.192857143 0.200000000 0.207142857 0.214285714
## [31] 0.221428571 0.228571429 0.235714286 0.242857143 0.250000000 0.257142857
## [37] 0.264285714 0.271428571 0.278571429 0.285714286 0.292857143 0.300000000
## [43] 0.307142857 0.314285714 0.321428571 0.328571429 0.335714286 0.342857143
## [49] 0.350000000 0.357142857 0.364285714 0.371428571 0.378571429 0.385714286
## [55] 0.392857143 0.400000000 0.407142857 0.414285714 0.421428571 0.428571429
## [61] 0.435714286 0.442857143 0.450000000 0.457142857 0.464285714 0.471428571
## [67] 0.478571429 0.485714286 0.492857143 0.500000000
##
## $db
## [1] 16.2878976 10.4415950 7.1554969 3.7705527 2.0578902 0.1234240
## [7] -0.5169944 -1.9700088 -2.8143084 -3.7168092 -4.6795792 -5.6470962
## [13] -6.1898830 -6.6860633 -7.2809263 -8.0526052 -8.5006366 -9.0594294
## [19] -9.0333225 -9.5760043 -10.5643890 -10.6765472 -10.8001562 -11.2714177
## [25] -11.5969942 -11.8298074 -12.2706154 -12.3273315 -12.2723218 -12.5014426
## [31] -13.0782008 -13.3890896 -13.6865067 -13.7782521 -13.8926049 -14.1044169
## [37] -14.5933890 -14.6062284 -14.9343244 -14.6908666 -14.9900725 -15.0286385
## [43] -15.0745157 -15.2621307 -15.8627498 -15.5180614 -15.5731139 -15.8900815
## [49] -15.8263756 -16.2333546 -16.7253183 -16.4583612 -16.5350865 -16.3455847
## [55] -16.4265913 -16.7418384 -16.6561977 -16.6060443 -17.0874251 -16.6036743
## [61] -16.3135863 -16.7180409 -16.8304366 -16.5422523 -16.6630950 -16.8029121
## [67] -16.9052827 -16.9961975 -16.7795028 -16.9091522
##
## $dbz
## [1] 11.6032654 11.0994077 10.2536556 9.0584317 7.5072577 5.6029201
## [7] 3.3763138 0.9242989 -1.5398967 -3.7018408 -5.3434863 -6.5237485
## [13] -7.4137624 -8.0795777 -8.5207761 -8.8075037 -9.0909353 -9.5062863
## [19] -10.1040273 -10.8405964 -11.6060119 -12.2840372 -12.8205859 -13.2360220
## [25] -13.5707120 -13.8378786 -14.0327139 -14.1710159 -14.3033983 -14.4891129
## [31] -14.7599543 -15.1038597 -15.4746944 -15.8208798 -16.1142639 -16.3563603
## [37] -16.5598162 -16.7289197 -16.8596253 -16.9549435 -17.0352369 -17.1305145
## [43] -17.2622796 -17.4312420 -17.6195370 -17.8045226 -17.9731102 -18.1252641
## [49] -18.2654994 -18.3927447 -18.4992085 -18.5787292 -18.6352288 -18.6818469
## [55] -18.7312105 -18.7860116 -18.8380398 -18.8760535 -18.8955282 -18.9019695
## [61] -18.9057367 -18.9141593 -18.9280759 -18.9446880 -18.9627140 -18.9842696
## [67] -19.0116116 -19.0422180 -19.0675217 -19.0774545plotts.sample.wge(as2$negative/us2$positive)## Warning in as2$negative/us2$positive: longer object length is not a multiple of
## shorter object length
## $autplt
## [1] 1.000000000 0.920600654 0.886730535 0.803011689 0.737135375
## [6] 0.671037046 0.605329951 0.544638239 0.486064603 0.424396460
## [11] 0.373448975 0.320146880 0.274513257 0.239068758 0.200964631
## [16] 0.169881086 0.133591348 0.107123767 0.082710921 0.062940720
## [21] 0.047152750 0.031579216 0.016492915 0.002446645 -0.005155609
## [26] -0.012435025
##
## $freq
## [1] 0.005494505 0.010989011 0.016483516 0.021978022 0.027472527 0.032967033
## [7] 0.038461538 0.043956044 0.049450549 0.054945055 0.060439560 0.065934066
## [13] 0.071428571 0.076923077 0.082417582 0.087912088 0.093406593 0.098901099
## [19] 0.104395604 0.109890110 0.115384615 0.120879121 0.126373626 0.131868132
## [25] 0.137362637 0.142857143 0.148351648 0.153846154 0.159340659 0.164835165
## [31] 0.170329670 0.175824176 0.181318681 0.186813187 0.192307692 0.197802198
## [37] 0.203296703 0.208791209 0.214285714 0.219780220 0.225274725 0.230769231
## [43] 0.236263736 0.241758242 0.247252747 0.252747253 0.258241758 0.263736264
## [49] 0.269230769 0.274725275 0.280219780 0.285714286 0.291208791 0.296703297
## [55] 0.302197802 0.307692308 0.313186813 0.318681319 0.324175824 0.329670330
## [61] 0.335164835 0.340659341 0.346153846 0.351648352 0.357142857 0.362637363
## [67] 0.368131868 0.373626374 0.379120879 0.384615385 0.390109890 0.395604396
## [73] 0.401098901 0.406593407 0.412087912 0.417582418 0.423076923 0.428571429
## [79] 0.434065934 0.439560440 0.445054945 0.450549451 0.456043956 0.461538462
## [85] 0.467032967 0.472527473 0.478021978 0.483516484 0.489010989 0.494505495
## [91] 0.500000000
##
## $db
## [1] 13.7596347 12.6170640 10.9202985 9.1159890 7.5096829 6.1937811
## [7] 5.0238297 3.7773446 2.4128435 0.9809050 -0.3825718 -1.4344221
## [13] -2.1026037 -2.5061577 -2.8357928 -3.3501795 -4.1160055 -5.0399034
## [19] -5.9834809 -6.6611745 -6.9052147 -6.8120137 -6.5646904 -6.4626854
## [25] -6.7349032 -7.3928020 -8.3777540 -9.4983035 -10.4914231 -11.2886558
## [31] -11.8625893 -12.2314873 -12.4987602 -12.6151821 -12.5568315 -12.4221080
## [37] -12.2672709 -12.2280247 -12.4192751 -12.7981588 -13.4070457 -14.2428783
## [43] -15.0948534 -15.7426842 -15.9702389 -15.7390035 -15.3434633 -14.9844708
## [49] -14.7657302 -14.7289775 -14.7027023 -14.5839269 -14.4593024 -14.4052325
## [55] -14.5493751 -14.9220179 -15.3550258 -15.7338895 -15.9037504 -15.6957188
## [61] -15.2396857 -14.7384552 -14.3191502 -14.0479295 -13.8292109 -13.5674287
## [67] -13.2212205 -12.6772192 -11.9302489 -11.1577653 -10.5402622 -10.2457280
## [73] -10.3100385 -10.6334680 -11.1023248 -11.5421142 -11.8059440 -11.8719254
## [79] -11.7371814 -11.4664057 -11.1724382 -10.8650390 -10.5501470 -10.2344037
## [85] -9.9232514 -9.7218241 -9.7083240 -9.8855850 -10.2198239 -10.5440572
## [91] -10.6705646
##
## $dbz
## [1] 10.7846516 10.5506095 10.1618961 9.6209367 8.9319557 8.1019108
## [7] 7.1418041 6.0683407 4.9056872 3.6866568 2.4520886 1.2469799
## [13] 0.1130318 -0.9199953 -1.8394046 -2.6479992 -3.3576252 -3.9827679
## [19] -4.5375615 -5.0359382 -5.4926732 -5.9235531 -6.3444004 -6.7696737
## [25] -7.2113160 -7.6779476 -8.1740402 -8.6987126 -9.2442160 -9.7948606
## [31] -10.3277975 -10.8171372 -11.2415203 -11.5923820 -11.8780739 -12.1208361
## [37] -12.3484507 -12.5853509 -12.8467412 -13.1364326 -13.4474256 -13.7641284
## [43] -14.0655910 -14.3295412 -14.5370102 -14.6768625 -14.7489241 -14.7643266
## [49] -14.7426513 -14.7069174 -14.6782599 -14.6717994 -14.6943135 -14.7436485
## [55] -14.8096004 -14.8760219 -14.9238520 -14.9345095 -14.8928293 -14.7888249
## [61] -14.6181431 -14.3817408 -14.0854754 -13.7398599 -13.3596796 -12.9630336
## [67] -12.5696948 -12.1991074 -11.8684809 -11.5912863 -11.3762033 -11.2264025
## [73] -11.1390419 -11.1049965 -11.1090659 -11.1311346 -11.1487905 -11.1414693
## [79] -11.0952229 -11.0062474 -10.8813877 -10.7353019 -10.5857176 -10.4488200
## [85] -10.3361232 -10.2531380 -10.1995344 -10.1703704 -10.1580434 -10.1546360
## [91] -10.1542296We can see that both the positive and negative cases are cumulative variables as well. However, looking at the ratio between the two, we can see that in Arizona the number of negative cases is vastly outweighing the positive cases.
Next lets look at the non cumulative data to see how it looks.
# This plot shows the increase in positive cases over time.
plotts.sample.wge(as2$positiveIncrease)
## $autplt
## [1] 1.0000000 0.8470508 0.8478437 0.8429955 0.8149689 0.8379732 0.7821325
## [8] 0.8185083 0.7573657 0.7408404 0.7030482 0.6896101 0.6698512 0.6416882
## [15] 0.6563625 0.5662050 0.5518993 0.5159686 0.4899257 0.4609225 0.4480741
## [22] 0.4340751 0.3323729 0.3404816 0.3067256 0.2951265
##
## $freq
## [1] 0.007142857 0.014285714 0.021428571 0.028571429 0.035714286 0.042857143
## [7] 0.050000000 0.057142857 0.064285714 0.071428571 0.078571429 0.085714286
## [13] 0.092857143 0.100000000 0.107142857 0.114285714 0.121428571 0.128571429
## [19] 0.135714286 0.142857143 0.150000000 0.157142857 0.164285714 0.171428571
## [25] 0.178571429 0.185714286 0.192857143 0.200000000 0.207142857 0.214285714
## [31] 0.221428571 0.228571429 0.235714286 0.242857143 0.250000000 0.257142857
## [37] 0.264285714 0.271428571 0.278571429 0.285714286 0.292857143 0.300000000
## [43] 0.307142857 0.314285714 0.321428571 0.328571429 0.335714286 0.342857143
## [49] 0.350000000 0.357142857 0.364285714 0.371428571 0.378571429 0.385714286
## [55] 0.392857143 0.400000000 0.407142857 0.414285714 0.421428571 0.428571429
## [61] 0.435714286 0.442857143 0.450000000 0.457142857 0.464285714 0.471428571
## [67] 0.478571429 0.485714286 0.492857143 0.500000000
##
## $db
## [1] 15.9206288 11.9774393 5.1441137 -0.8082256 -5.2527223 -4.0931875
## [7] -5.7085142 -10.2790733 -4.2037479 -4.5663927 -7.0051850 -15.4673458
## [13] -9.4187946 -7.4354190 -7.4054209 -10.1420896 -13.5714321 -17.6425756
## [19] -12.2313254 -3.8192633 -3.9299655 -4.6767915 -8.9486672 -12.0583288
## [25] -12.9107780 -21.8587070 -12.4310039 -7.1817034 -5.1049432 -8.3391835
## [31] -14.2919696 -14.1516361 -9.7705523 -5.7673327 -6.0308224 -9.4731382
## [37] -21.4545114 -18.1866036 -8.4441442 -3.1471867 -3.7526140 -8.8688047
## [43] -15.4381695 -12.9424401 -20.5755347 -14.3359730 -13.7478544 -7.5465057
## [49] -5.3644482 -6.1643723 -12.3712356 -13.4826921 -8.7681009 -5.2323980
## [55] -5.0440218 -8.8583051 -16.0575587 -8.5162585 -4.2059821 -1.6660321
## [61] -1.7721567 -6.7568759 -15.1738396 -11.5667492 -12.0825769 -17.4122063
## [67] -14.9982751 -14.7491544 -17.3452283 -25.7702284
##
## $dbz
## [1] 11.2486368 10.7378757 9.8774417 8.6536857 7.0483872 5.0421055
## [7] 2.6252550 -0.1689095 -3.1795772 -5.9827190 -7.9754489 -8.9683920
## [13] -9.2844931 -9.1763152 -8.7510138 -8.1569966 -7.5736868 -7.1248268
## [19] -6.8613385 -6.7895653 -6.8991493 -7.1753902 -7.5957207 -8.1152556
## [25] -8.6519878 -9.0916440 -9.3326781 -9.3507981 -9.2130240 -9.0216391
## [31] -8.8525672 -8.7329379 -8.6491605 -8.5671656 -8.4555462 -8.3048241
## [37] -8.1351888 -7.9900387 -7.9211618 -7.9737653 -8.1746435 -8.5214572
## [43] -8.9703685 -9.4258540 -9.7520721 -9.8309888 -9.6415220 -9.2713156
## [49] -8.8444242 -8.4504357 -8.1220857 -7.8443067 -7.5746104 -7.2689841
## [55] -6.9079839 -6.5107370 -6.1286563 -5.8264375 -5.6650430 -5.6937727
## [61] -5.9495701 -6.4591354 -7.2405428 -8.3021694 -9.6365303 -11.2050965
## [67] -12.9087869 -14.5447781 -15.7828924 -16.2533971# This shows the same for negative cases over time.
plotts.sample.wge(as2$negativeIncrease)
## $autplt
## [1] 1.0000000 0.6939751 0.7510290 0.7357039 0.6767746 0.7549883 0.6873635
## [8] 0.6870708 0.6326626 0.6486209 0.6282344 0.5800424 0.5548354 0.5773257
## [15] 0.6125914 0.5273803 0.5488208 0.4454643 0.4667526 0.5163171 0.4349937
## [22] 0.4641381 0.3766452 0.4018404 0.3765581 0.3680332
##
## $freq
## [1] 0.007142857 0.014285714 0.021428571 0.028571429 0.035714286 0.042857143
## [7] 0.050000000 0.057142857 0.064285714 0.071428571 0.078571429 0.085714286
## [13] 0.092857143 0.100000000 0.107142857 0.114285714 0.121428571 0.128571429
## [19] 0.135714286 0.142857143 0.150000000 0.157142857 0.164285714 0.171428571
## [25] 0.178571429 0.185714286 0.192857143 0.200000000 0.207142857 0.214285714
## [31] 0.221428571 0.228571429 0.235714286 0.242857143 0.250000000 0.257142857
## [37] 0.264285714 0.271428571 0.278571429 0.285714286 0.292857143 0.300000000
## [43] 0.307142857 0.314285714 0.321428571 0.328571429 0.335714286 0.342857143
## [49] 0.350000000 0.357142857 0.364285714 0.371428571 0.378571429 0.385714286
## [55] 0.392857143 0.400000000 0.407142857 0.414285714 0.421428571 0.428571429
## [61] 0.435714286 0.442857143 0.450000000 0.457142857 0.464285714 0.471428571
## [67] 0.478571429 0.485714286 0.492857143 0.500000000
##
## $db
## [1] 15.7817498 8.0390858 5.9061951 3.6757892 -4.2416109 -0.4532261
## [7] -2.3437355 -8.4204904 -0.7359364 -3.7533347 -10.6901767 -14.0850224
## [13] -27.4387203 -11.7175911 -11.5514070 -5.7490633 -17.7695717 -7.0309068
## [19] -5.0063359 -2.4389114 -3.3575708 -4.7137641 -10.3003869 -14.2765601
## [25] -11.2684307 -20.6785433 -10.0149714 -0.8116063 -0.7359700 -2.9762084
## [31] -7.8928165 -17.0204960 -18.8370627 -15.7137330 -12.7929001 -14.7856460
## [37] -10.0801750 -23.4393600 -8.8912824 -8.9877847 -4.9872271 -4.4838716
## [43] -7.8446433 -2.4601378 -7.0097202 -3.3559861 -7.6389939 -24.1230678
## [49] -9.4373421 -4.2895550 -1.0196580 -12.1087946 -3.8761438 1.5808958
## [55] -14.9896665 -5.1932899 -17.9058113 -11.1436147 -0.6133112 1.4821512
## [61] -0.6967830 -9.3275365 -16.5475593 -3.6107026 -5.8379556 -5.0318682
## [67] -14.0169902 -11.5481727 -3.6728605 -14.2680347
##
## $dbz
## [1] 10.6568176 10.1446545 9.2845084 8.0680435 6.4880152 4.5472676
## [7] 2.2798615 -0.2062824 -2.6750022 -4.7905686 -6.3416527 -7.3834934
## [13] -8.0030691 -8.1502175 -7.8223089 -7.2179217 -6.5956804 -6.1200606
## [19] -5.8506590 -5.7829142 -5.8758627 -6.0606871 -6.2435232 -6.3236380
## [25] -6.2388860 -6.0098179 -5.7330243 -5.5295000 -5.4996361 -5.7075191
## [31] -6.1814847 -6.9132519 -7.8441651 -8.8341993 -9.6375568 -9.9744336
## [37] -9.7395549 -9.0946332 -8.2934237 -7.5186194 -6.8627224 -6.3603344
## [43] -6.0114482 -5.7901371 -5.6473419 -5.5176787 -5.3378450 -5.0734889
## [49] -4.7359697 -4.3727596 -4.0395047 -3.7750465 -3.5902994 -3.4698155
## [55] -3.3827152 -3.3001635 -3.2130590 -3.1400612 -3.1210444 -3.2016762
## [61] -3.4186463 -3.7898719 -4.3083179 -4.9372935 -5.6089912 -6.2340376
## [67] -6.7296420 -7.0555684 -7.2250938 -7.2753982The positive rate shows a bit of a lag to the US data, which shows a sharp increase in April, to a decrease in May-June to another increase in July. The Negatives however, gradually increase as a relationship to the increase in testing that went on in the state.
The totalIncrease created by summing the positive and negative increase should be equivalent to the number of tests administered without considering pending cases.
totalIncrease = as2$positiveIncrease + as2$negativeIncreaseThis variable (totalTestResultsIncresase) should closely match the totalIncrease
plotts.sample.wge(totalIncrease)
## $autplt
## [1] 1.0000000 0.7487270 0.7936271 0.7748428 0.7284560 0.7974090 0.7265687
## [8] 0.7369464 0.6778396 0.6877675 0.6619148 0.6286897 0.6023545 0.6119466
## [15] 0.6463648 0.5537321 0.5680936 0.4809054 0.4931735 0.5235945 0.4581411
## [22] 0.4780621 0.3847051 0.4086774 0.3760781 0.3691591
##
## $freq
## [1] 0.007142857 0.014285714 0.021428571 0.028571429 0.035714286 0.042857143
## [7] 0.050000000 0.057142857 0.064285714 0.071428571 0.078571429 0.085714286
## [13] 0.092857143 0.100000000 0.107142857 0.114285714 0.121428571 0.128571429
## [19] 0.135714286 0.142857143 0.150000000 0.157142857 0.164285714 0.171428571
## [25] 0.178571429 0.185714286 0.192857143 0.200000000 0.207142857 0.214285714
## [31] 0.221428571 0.228571429 0.235714286 0.242857143 0.250000000 0.257142857
## [37] 0.264285714 0.271428571 0.278571429 0.285714286 0.292857143 0.300000000
## [43] 0.307142857 0.314285714 0.321428571 0.328571429 0.335714286 0.342857143
## [49] 0.350000000 0.357142857 0.364285714 0.371428571 0.378571429 0.385714286
## [55] 0.392857143 0.400000000 0.407142857 0.414285714 0.421428571 0.428571429
## [61] 0.435714286 0.442857143 0.450000000 0.457142857 0.464285714 0.471428571
## [67] 0.478571429 0.485714286 0.492857143 0.500000000
##
## $db
## [1] 15.9085677 9.4118568 5.9486713 2.7457293 -4.2610083 -1.2110685
## [7] -3.0204479 -8.8141476 -1.4274239 -3.7512935 -9.4439698 -18.4133153
## [13] -18.9540627 -10.5833759 -10.1251657 -8.3377816 -23.9923378 -9.9995727
## [19] -6.5183938 -3.0152818 -3.3212456 -4.5508449 -9.9698371 -13.6594802
## [25] -11.6664202 -21.3193922 -10.3281565 -1.8622303 -1.4289392 -3.8980239
## [31] -9.2149095 -16.0062412 -15.8993205 -14.3144850 -11.9003456 -13.4055037
## [37] -13.1177672 -23.2562918 -8.5553132 -8.4640344 -4.8670927 -5.5500905
## [43] -9.0624129 -4.5491781 -9.8608665 -4.9033679 -10.0960996 -16.3989455
## [49] -8.0559367 -4.7832459 -2.5171243 -12.2835385 -4.7379287 0.5513428
## [55] -11.0788915 -6.5678097 -21.1949817 -12.5829978 -1.1921107 0.9516414
## [61] -0.7370283 -8.4090350 -16.0882235 -4.7759968 -7.0081127 -6.8908509
## [67] -14.0217439 -13.3254458 -5.4290272 -15.7760063
##
## $dbz
## [1] 10.9057384 10.3909614 9.5253468 8.2984014 6.6986362 4.7207277
## [7] 2.3837375 -0.2276968 -2.8960927 -5.2570111 -7.0161916 -8.1946468
## [13] -8.9088059 -9.1058428 -8.7545096 -8.0686846 -7.3473432 -6.7770556
## [19] -6.4217599 -6.2770244 -6.3062070 -6.4520908 -6.6370870 -6.7706230
## [25] -6.7801291 -6.6551419 -6.4615968 -6.3074452 -6.2952609 -6.4961734
## [31] -6.9436033 -7.6301247 -8.4944879 -9.3964731 -10.1077693 -10.3946021
## [37] -10.1926561 -9.6560387 -9.0041627 -8.3932112 -7.9004922 -7.5493611
## [43] -7.3276270 -7.1946874 -7.0864372 -6.9293725 -6.6687864 -6.2971271
## [49] -5.8571194 -5.4153486 -5.0295630 -4.7295079 -4.5130558 -4.3527990
## [55] -4.2107625 -4.0587457 -3.8959345 -3.7528182 -3.6795464 -3.7284158
## [61] -3.9404330 -4.3380703 -4.9208380 -5.6600575 -6.4927087 -7.3212047
## [67] -8.0329803 -8.5438280 -8.8323493 -8.9227701The totalincrease variable showed to have a similar pattern as the positive cases for the arizona data. The totalIncrease variable shows to match exactly with the totalTestResultsIncrease variable, so now we can have confidence on how that variable was calculated and use it moving forward.
Here is a graph of increase count in total tests.
plotts.sample.wge(as2$totalTestResultsIncrease)
## $autplt
## [1] 1.0000000 0.7487270 0.7936271 0.7748428 0.7284560 0.7974090 0.7265687
## [8] 0.7369464 0.6778396 0.6877675 0.6619148 0.6286897 0.6023545 0.6119466
## [15] 0.6463648 0.5537321 0.5680936 0.4809054 0.4931735 0.5235945 0.4581411
## [22] 0.4780621 0.3847051 0.4086774 0.3760781 0.3691591
##
## $freq
## [1] 0.007142857 0.014285714 0.021428571 0.028571429 0.035714286 0.042857143
## [7] 0.050000000 0.057142857 0.064285714 0.071428571 0.078571429 0.085714286
## [13] 0.092857143 0.100000000 0.107142857 0.114285714 0.121428571 0.128571429
## [19] 0.135714286 0.142857143 0.150000000 0.157142857 0.164285714 0.171428571
## [25] 0.178571429 0.185714286 0.192857143 0.200000000 0.207142857 0.214285714
## [31] 0.221428571 0.228571429 0.235714286 0.242857143 0.250000000 0.257142857
## [37] 0.264285714 0.271428571 0.278571429 0.285714286 0.292857143 0.300000000
## [43] 0.307142857 0.314285714 0.321428571 0.328571429 0.335714286 0.342857143
## [49] 0.350000000 0.357142857 0.364285714 0.371428571 0.378571429 0.385714286
## [55] 0.392857143 0.400000000 0.407142857 0.414285714 0.421428571 0.428571429
## [61] 0.435714286 0.442857143 0.450000000 0.457142857 0.464285714 0.471428571
## [67] 0.478571429 0.485714286 0.492857143 0.500000000
##
## $db
## [1] 15.9085677 9.4118568 5.9486713 2.7457293 -4.2610083 -1.2110685
## [7] -3.0204479 -8.8141476 -1.4274239 -3.7512935 -9.4439698 -18.4133153
## [13] -18.9540627 -10.5833759 -10.1251657 -8.3377816 -23.9923378 -9.9995727
## [19] -6.5183938 -3.0152818 -3.3212456 -4.5508449 -9.9698371 -13.6594802
## [25] -11.6664202 -21.3193922 -10.3281565 -1.8622303 -1.4289392 -3.8980239
## [31] -9.2149095 -16.0062412 -15.8993205 -14.3144850 -11.9003456 -13.4055037
## [37] -13.1177672 -23.2562918 -8.5553132 -8.4640344 -4.8670927 -5.5500905
## [43] -9.0624129 -4.5491781 -9.8608665 -4.9033679 -10.0960996 -16.3989455
## [49] -8.0559367 -4.7832459 -2.5171243 -12.2835385 -4.7379287 0.5513428
## [55] -11.0788915 -6.5678097 -21.1949817 -12.5829978 -1.1921107 0.9516414
## [61] -0.7370283 -8.4090350 -16.0882235 -4.7759968 -7.0081127 -6.8908509
## [67] -14.0217439 -13.3254458 -5.4290272 -15.7760063
##
## $dbz
## [1] 10.9057384 10.3909614 9.5253468 8.2984014 6.6986362 4.7207277
## [7] 2.3837375 -0.2276968 -2.8960927 -5.2570111 -7.0161916 -8.1946468
## [13] -8.9088059 -9.1058428 -8.7545096 -8.0686846 -7.3473432 -6.7770556
## [19] -6.4217599 -6.2770244 -6.3062070 -6.4520908 -6.6370870 -6.7706230
## [25] -6.7801291 -6.6551419 -6.4615968 -6.3074452 -6.2952609 -6.4961734
## [31] -6.9436033 -7.6301247 -8.4944879 -9.3964731 -10.1077693 -10.3946021
## [37] -10.1926561 -9.6560387 -9.0041627 -8.3932112 -7.9004922 -7.5493611
## [43] -7.3276270 -7.1946874 -7.0864372 -6.9293725 -6.6687864 -6.2971271
## [49] -5.8571194 -5.4153486 -5.0295630 -4.7295079 -4.5130558 -4.3527990
## [55] -4.2107625 -4.0587457 -3.8959345 -3.7528182 -3.6795464 -3.7284158
## [61] -3.9404330 -4.3380703 -4.9208380 -5.6600575 -6.4927087 -7.3212047
## [67] -8.0329803 -8.5438280 -8.8323493 -8.9227701Now we are divding the total test by cumulative positive cases.
pp = totalIncrease / as2$positive
plotts.sample.wge(pp[18:140])
## $autplt
## [1] 1.0000000000 0.1009224279 0.4908693807 0.2058653775 0.1681307962
## [6] 0.1216669929 0.1181007607 0.0844962626 0.1538612397 -0.0019239054
## [11] 0.1425198193 0.0290542918 0.0131829181 0.0635085534 -0.0015806895
## [16] 0.0359868231 0.0078682229 0.0033079042 -0.0005962703 0.0033669252
## [21] 0.0138206571 0.0060047184 -0.0028636199 -0.0071876085 -0.0065889789
## [26] -0.0151612007
##
## $freq
## [1] 0.008130081 0.016260163 0.024390244 0.032520325 0.040650407 0.048780488
## [7] 0.056910569 0.065040650 0.073170732 0.081300813 0.089430894 0.097560976
## [13] 0.105691057 0.113821138 0.121951220 0.130081301 0.138211382 0.146341463
## [19] 0.154471545 0.162601626 0.170731707 0.178861789 0.186991870 0.195121951
## [25] 0.203252033 0.211382114 0.219512195 0.227642276 0.235772358 0.243902439
## [31] 0.252032520 0.260162602 0.268292683 0.276422764 0.284552846 0.292682927
## [37] 0.300813008 0.308943089 0.317073171 0.325203252 0.333333333 0.341463415
## [43] 0.349593496 0.357723577 0.365853659 0.373983740 0.382113821 0.390243902
## [49] 0.398373984 0.406504065 0.414634146 0.422764228 0.430894309 0.439024390
## [55] 0.447154472 0.455284553 0.463414634 0.471544715 0.479674797 0.487804878
## [61] 0.495934959
##
## $db
## [1] 6.07542576 6.75118438 6.43629420 4.21183897 3.29674558 2.92132968
## [7] 1.18921208 0.81268091 0.05555239 -0.17019597 -0.56362484 -0.88196065
## [13] -1.15131345 -1.00543710 -2.86396700 -2.03967242 -3.57678411 -1.36351058
## [19] -4.67423771 -6.93691543 -6.43689356 -6.62578927 -5.97092540 -5.58768881
## [25] -6.94898216 -6.71084915 -7.60233605 -8.99137676 -8.91880097 -8.44121233
## [31] -6.84409244 -5.09089375 -5.26114173 -6.14751680 -2.79127484 -2.33427197
## [37] -2.24759992 -2.50530961 -3.72533735 -2.48279011 -5.13164293 -3.44050101
## [43] -3.87388511 -1.71427235 0.09528964 0.96017613 1.37509959 2.03164756
## [49] 2.17270193 0.94100261 0.88468732 -0.98396938 -0.25468087 -0.82899025
## [55] -0.32253609 1.07820186 1.40733127 2.50077101 2.27338176 2.77475806
## [61] 2.97353977
##
## $dbz
## [1] 5.38054316 5.13207247 4.73213490 4.20326665 3.57865772 2.90141510
## [7] 2.22010919 1.57952790 1.00907677 0.51509199 0.08224981 -0.31737648
## [13] -0.71406519 -1.13395740 -1.59675779 -2.11533452 -2.69488931 -3.33069107
## [19] -4.00506554 -4.68612335 -5.33194462 -5.90231963 -6.37342983 -6.74371980
## [25] -7.02293562 -7.21089368 -7.28268479 -7.19478938 -6.91372812 -6.44761815
## [31] -5.85202749 -5.20657243 -4.58647258 -4.04745402 -3.62299775 -3.32566977
## [37] -3.14689579 -3.05380054 -2.98605184 -2.86115046 -2.59882238 -2.16090472
## [43] -1.57625617 -0.92466012 -0.29694206 0.23470814 0.62620482 0.85998370
## [49] 0.94179437 0.89954951 0.78323870 0.66126474 0.60679073 0.67309998
## [55] 0.87102500 1.16738713 1.50564166 1.83029054 2.10021358 2.29038585
## [61] 2.38800552Now we are creating the variable posper, which stands for positive percentage, given as decimal percentages.
posperas = (as2$positiveIncrease[50:140]/as2$totalTestResultsIncrease[50:140])
plotts.sample.wge(posperas)
## $autplt
## [1] 1.00000000 0.71792171 0.67922925 0.76957696 0.66643001 0.62690915
## [7] 0.68724945 0.67077652 0.64696508 0.63459787 0.65448167 0.55874413
## [13] 0.47073896 0.50573985 0.45471971 0.41128371 0.43264065 0.36214773
## [19] 0.32279291 0.33135597 0.23047029 0.22112612 0.16296073 0.10908299
## [25] 0.09953088 0.08399247
##
## $freq
## [1] 0.01098901 0.02197802 0.03296703 0.04395604 0.05494505 0.06593407
## [7] 0.07692308 0.08791209 0.09890110 0.10989011 0.12087912 0.13186813
## [13] 0.14285714 0.15384615 0.16483516 0.17582418 0.18681319 0.19780220
## [19] 0.20879121 0.21978022 0.23076923 0.24175824 0.25274725 0.26373626
## [25] 0.27472527 0.28571429 0.29670330 0.30769231 0.31868132 0.32967033
## [31] 0.34065934 0.35164835 0.36263736 0.37362637 0.38461538 0.39560440
## [37] 0.40659341 0.41758242 0.42857143 0.43956044 0.45054945 0.46153846
## [43] 0.47252747 0.48351648 0.49450549
##
## $db
## [1] 15.2350205 -7.1424843 -28.4923951 -19.2285775 -9.6015452 -8.0580807
## [7] -6.4997626 -5.8492410 -1.0937259 0.4416191 -10.7729372 -5.7528722
## [13] -5.4577631 -25.7025087 -9.4393295 -15.8528829 -17.3372605 -2.1238730
## [19] -8.4545922 -10.9650402 -10.5131139 -11.0845339 -13.1491509 -7.8652841
## [25] -4.6057396 1.0031854 -7.8280142 0.4581398 -1.6291984 -5.4699216
## [31] -6.4217203 -8.0863315 -4.3654074 -1.4972342 -3.7132038 -5.1565554
## [37] -10.3589614 -6.2548481 -6.2141164 -14.4889082 -18.7173614 -14.2235246
## [43] -7.7095331 -4.9035540 -11.0082189
##
## $dbz
## [1] 9.7019918 8.8667680 7.4523014 5.4332491 2.8145574 -0.2269771
## [7] -2.9735427 -4.2889653 -4.3394497 -4.1856901 -4.2985220 -4.7363983
## [13] -5.4514855 -6.3383292 -7.1857873 -7.7129813 -7.8026298 -7.6306599
## [19] -7.4350128 -7.2852997 -7.0626970 -6.5851960 -5.7991817 -4.8456733
## [25] -3.9345003 -3.2223098 -2.7874386 -2.6456960 -2.7632524 -3.0630783
## [31] -3.4406820 -3.8054890 -4.1333109 -4.4792800 -4.9325737 -5.5554368
## [37] -6.3418942 -7.1973298 -7.9494217 -8.4269822 -8.5804295 -8.5057221
## [43] -8.3443401 -8.1979205 -8.1161032Looking at the Arizona positive percentage indicates that the rate of positive tests and the likelihood of a test to be positive increased over time through June and looks like it is beginning to dip in July.
AZ Forecast
Stationary Forecast
aic5.wge(posperas, p=0:5, q=0:2) # arma(3,0)## ---------WORKING... PLEASE WAIT...
##
##
## Five Smallest Values of aic## p q aic
## 10 3 0 -5.956097
## 11 3 1 -5.940761
## 14 4 1 -5.939006
## 13 4 0 -5.937297
## 12 3 2 -5.935570aic5.wge(posperas, p=0:5, q=0:2, type='bic') # arma(1,1)## ---------WORKING... PLEASE WAIT...
##
##
## Five Smallest Values of bic## p q bic
## 5 1 1 -5.846965
## 10 3 0 -5.845730
## 11 3 1 -5.802802
## 6 1 2 -5.802276
## 8 2 1 -5.799875The BIC model selects the ARMA (1,1) model which we will use moving forward.
estPosper22 = est.arma.wge(posperas, p=1, q=1)##
## Coefficients of Original polynomial:
## 0.9805
##
## Factor Roots Abs Recip System Freq
## 1-0.9805B 1.0198 0.9805 0.0000
##
## mean(posperas)## [1] 0.1647793estPosper22$phi## [1] 0.9805499estPosper22$theta## [1] 0.7312676After running our estimated model based on an ARMA(1,1) model, we end up with the model below.
(1 - 0.98B)(X - 0.16) = (1 - 0.73a_t)
factor.wge(phi=estPosper22$phi)##
## Coefficients of Original polynomial:
## 0.9805
##
## Factor Roots Abs Recip System Freq
## 1-0.9805B 1.0198 0.9805 0.0000
##
## Forcasting 10 ahead for the ARMA model. Both of these will eventually reach the mean of 0.16
fore.arma.wge(posperas,phi=estPosper22$phi, theta=estPosper22$theta, n.ahead=10,plot=TRUE)
## $f
## [1] 0.2355517 0.2341752 0.2328254 0.2315019 0.2302041 0.2289316 0.2276839
## [8] 0.2264604 0.2252606 0.2240843
##
## $ll
## [1] 0.1372237 0.1328380 0.1286770 0.1247204 0.1209507 0.1173531 0.1139147
## [8] 0.1106239 0.1074708 0.1044464
##
## $ul
## [1] 0.3338797 0.3355123 0.3369738 0.3382834 0.3394576 0.3405101 0.3414530
## [8] 0.3422968 0.3430505 0.3437222
##
## $resid
## [1] 0.000000000 0.003940524 -0.008510589 -0.012611691 -0.034054685
## [6] -0.023147882 0.131519510 0.023522190 -0.017121450 -0.061304751
## [11] -0.018651346 -0.054860352 -0.062518842 0.013267722 -0.003756197
## [16] -0.068649269 0.020763138 -0.016830704 -0.071532037 -0.046773591
## [21] -0.023636393 0.014428823 -0.004967516 -0.016363160 0.018000853
## [26] -0.023143488 -0.031016305 0.003304884 0.008707622 -0.013327293
## [31] -0.012188119 0.019388032 0.009394469 -0.034593999 -0.005790202
## [36] 0.005285481 -0.008835845 0.018192305 0.019530240 -0.012539478
## [41] -0.028924078 0.027056279 0.050281627 -0.011044014 0.003132813
## [46] 0.033167942 0.024195374 -0.004731925 -0.006005050 0.114671037
## [51] -0.012322674 -0.026807467 0.128416693 -0.050286764 -0.024691867
## [56] 0.104397517 -0.031199884 0.008970941 0.080931691 0.039004448
## [61] -0.028832209 0.018587418 0.130962593 -0.028804303 -0.023652784
## [66] 0.004957484 0.014398351 0.031771121 0.179675310 -0.049486629
## [71] 0.024180032 0.033541805 -0.016654274 -0.020972652 -0.017984002
## [76] -0.001093209 0.085573980 0.019241479 0.067972793 -0.023892841
## [81] -0.071393797 -0.048718604 -0.026057624 -0.016083561 0.002682210
## [86] 0.015911598 0.017244178 0.155667034 -0.103563674 -0.037647837
## [91] 0.014101782
##
## $wnv
## [1] 0.002516764
##
## $se
## [1] 0.05016736 0.05170261 0.05313691 0.05448037 0.05574154 0.05692780
## [7] 0.05804550 0.05910021 0.06009685 0.06103975
##
## $psi
## [1] 0.2492823 0.2444337 0.2396795 0.2350177 0.2304466 0.2259644 0.2215693
## [8] 0.2172598 0.2130341 0.2088905Looking at the forecast, we can see that with the AR1 model, our forecasts are regressing back to the mean.
Now let’s check 20 steps back.
for22 = fore.arma.wge(posperas,phi=estPosper22$phi, theta=estPosper22$theta, n.ahead=20,lastn=TRUE, plot=TRUE)
ASEa22 = mean((posperas[(90-20+1):90] - for22$f)^2)
ASEa22## [1] 0.002586599The ASE score for the stationary Arizona model was 0.0026 which is relatively low.
trainingSize = 70
horizon = 10
ASEHolder = numeric()
for( i in 1:(81-(trainingSize + horizon) + 1))
{
forecasts = fore.aruma.wge(posperas[i:(i+(trainingSize-1))],phi = estPosper22$phi, theta = estPosper22$theta, s = 0, d = 0,n.ahead = horizon)
ASE = mean((posperas[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
ASEHolder[i] = ASE
}
ASEHolder## [1] 0.002312777 0.002053826hist(ASEHolder)
WindowedASE = mean(ASEHolder)
summary(ASEHolder)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.002054 0.002119 0.002183 0.002183 0.002248 0.002313WindowedASE## [1] 0.002183301Running our rolling window average, our ASE holder is hovering around 0.0020 and 0.0023 within the 10 day lookback window.
Now lets try adding a nonstationary process.
p.ns = artrans.wge(posperas,phi.tr=(1))
plotts.sample.wge(p.ns)
## $autplt
## [1] 1.00000000 -0.43086310 -0.22552536 0.29867778 -0.08499842 -0.17995361
## [7] 0.15049317 0.01184274 -0.02250657 -0.05953430 0.21101065 -0.02951573
## [13] -0.24193439 0.17017731 -0.02571653 -0.10248400 0.17232707 -0.05959130
## [19] -0.09313106 0.18765138 -0.16252150 0.10312189 -0.05161131 -0.04655373
## [25] 0.01664135 0.08523630
##
## $freq
## [1] 0.01111111 0.02222222 0.03333333 0.04444444 0.05555556 0.06666667
## [7] 0.07777778 0.08888889 0.10000000 0.11111111 0.12222222 0.13333333
## [13] 0.14444444 0.15555556 0.16666667 0.17777778 0.18888889 0.20000000
## [19] 0.21111111 0.22222222 0.23333333 0.24444444 0.25555556 0.26666667
## [25] 0.27777778 0.28888889 0.30000000 0.31111111 0.32222222 0.33333333
## [31] 0.34444444 0.35555556 0.36666667 0.37777778 0.38888889 0.40000000
## [37] 0.41111111 0.42222222 0.43333333 0.44444444 0.45555556 0.46666667
## [43] 0.47777778 0.48888889 0.50000000
##
## $db
## [1] -8.4550001 -21.4574773 -16.0084325 -13.1759750 -23.7474942 -29.4410977
## [7] -19.5505314 -15.4892002 -5.0418335 -0.3483469 -6.3899708 -8.0990085
## [13] -3.8452171 -31.1346481 -5.7605806 -8.8502304 -17.5928154 2.6937865
## [19] 0.6314160 -4.9144582 -4.9698736 -1.3607796 -5.5857595 -17.3435722
## [25] -3.6474248 7.1832716 -1.0299279 8.3434123 5.1064913 -0.1046000
## [31] 0.6133391 -5.0262075 7.0810212 2.9894198 3.5761768 -4.9008519
## [37] -11.7551905 0.6801849 2.0539758 -1.7021376 -1.8911150 3.3450612
## [43] 2.5698362 -4.3507638 -0.5027917
##
## $dbz
## [1] -12.92660167 -13.23098552 -13.49574294 -13.26227875 -12.15736854
## [6] -10.41879417 -8.59402100 -7.03402529 -5.86919632 -5.11864359
## [11] -4.74393802 -4.65490049 -4.69795053 -4.67176459 -4.41975313
## [16] -3.94337331 -3.38040804 -2.86754337 -2.43997666 -2.00858315
## [21] -1.40830783 -0.52999046 0.56373901 1.68753280 2.66527566
## [26] 3.39305197 3.83546622 4.00820956 3.96640374 3.79295637
## [31] 3.57396579 3.36008924 3.13979120 2.85325293 2.43750151
## [36] 1.87147326 1.20674651 0.57646089 0.15190061 0.03238569
## [41] 0.16810080 0.41857753 0.65856023 0.81917912 0.87455983We accounted for the wandering behavior with our trend forecast. Which was a low pass filter.
acf(p.ns)
aic5.wge(p.ns, p=0:10, q=0:2) # arima(10,1,2)## ---------WORKING... PLEASE WAIT...
##
##
## Error in aic calculation at 5 2
## Five Smallest Values of aic## p q aic
## 33 10 2 -6.167285
## 32 10 1 -6.109982
## 21 6 2 -6.046182
## 31 10 0 -5.979502
## 26 8 1 -5.960030aic5.wge(p.ns, p=0:10, q=0:2, type='bic') #arima(2,1,0)## ---------WORKING... PLEASE WAIT...
##
##
## Error in aic calculation at 5 2
## Five Smallest Values of bic## p q bic
## 2 0 1 -5.892024
## 7 2 0 -5.864710
## 3 0 2 -5.845232
## 5 1 1 -5.843990
## 8 2 1 -5.841884Looking at the ACF of the differenced data, we can see that a couple of the residuals fall out of the 95% confidence intervals, which is not great. However, we can perform the Ljung Box tex test to see how it does.
ljung.wge(p.ns)## Obs -0.4308631 -0.2255254 0.2986778 -0.08499842 -0.1799536 0.1504932 0.01184274 -0.02250657 -0.0595343 0.2110107 -0.02951573 -0.2419344 0.1701773 -0.02571653 -0.102484 0.1723271 -0.0595913 -0.09313106 0.1876514 -0.1625215 0.1031219 -0.05161131 -0.04655373 0.01664135## $test
## [1] "Ljung-Box test"
##
## $K
## [1] 24
##
## $chi.square
## [1] 66.17265
##
## $df
## [1] 24
##
## $pval
## [1] 8.178399e-06ljung.wge(p.ns, K = 48)## Obs -0.4308631 -0.2255254 0.2986778 -0.08499842 -0.1799536 0.1504932 0.01184274 -0.02250657 -0.0595343 0.2110107 -0.02951573 -0.2419344 0.1701773 -0.02571653 -0.102484 0.1723271 -0.0595913 -0.09313106 0.1876514 -0.1625215 0.1031219 -0.05161131 -0.04655373 0.01664135 0.0852363 -0.04378138 -0.04856567 0.02753952 0.03993 -0.04546079 -0.02666046 0.08786457 -0.06172621 -0.1039269 0.1652229 -0.01748764 -0.1377251 0.130014 -0.01323147 -0.09400176 0.05898378 0.03096454 -0.03726662 -0.04811551 0.03908592 0.08999835 -0.1292868 0.00538244## $test
## [1] "Ljung-Box test"
##
## $K
## [1] 48
##
## $chi.square
## [1] 89.35886
##
## $df
## [1] 48
##
## $pval
## [1] 0.0002704864Running the Ljung box test, we can see that both p values are below 0.05, which fails to reject the null hypothesis of no white noise.
That said, we will proceed to see how our model performs.
Running our AIC, we can see that the preferred ARIMA model was (10,1,2) while our BIC preferred (2,1,0), which we will move forward with.
estp.ns = est.arma.wge(p.ns,p=2,q=0)##
## Coefficients of Original polynomial:
## -0.6509 -0.5049
##
## Factor Roots Abs Recip System Freq
## 1+0.6509B+0.5049B^2 -0.6446+-1.2511i 0.7105 0.3257
##
## The preferred estimated model is (1 + 0.65B + 0.50B^2)(X - B)(1 - B) = a_t
Using this model, we will forecast a nonstationary ARIMA model.
fore.arma.wge(p.ns,phi=estp.ns$phi, n.ahead=10,plot=TRUE,limits=FALSE)
## $f
## [1] -0.0432687722 0.0100449375 0.0177961903 -0.0141654941 0.0027247849
## [6] 0.0078674824 -0.0040072454 0.0011255466 0.0037798299 -0.0005392136
##
## $ll
## [1] -0.14121368 -0.10682029 -0.09933939 -0.13711574 -0.12189068 -0.11687590
## [7] -0.12952566 -0.12454657 -0.12192778 -0.12634948
##
## $ul
## [1] 0.05467613 0.12691017 0.13493177 0.10878476 0.12734024 0.13261087
## [7] 0.12151117 0.12679767 0.12948744 0.12527105
##
## $resid
## [1] 0.0000000000 0.0000000000 -0.0133374578 -0.0360543871 -0.0182928827
## [6] 0.1366998386 0.0217957512 -0.0101981230 -0.1111873731 -0.0241708442
## [11] -0.0496575141 -0.0359437382 0.0243155160 0.0138258461 -0.0449991809
## [16] 0.0219839493 -0.0188621055 -0.0436819520 -0.0478224463 -0.0134944459
## [21] 0.0438359322 0.0120068963 -0.0055540806 0.0151715953 -0.0222464554
## [26] -0.0212279968 0.0002146923 0.0177690971 -0.0011958820 -0.0108259214
## [31] 0.0181895260 0.0134239289 -0.0295791086 -0.0086593446 0.0026399186
## [36] 0.0046011229 0.0223121514 0.0164982552 -0.0100538982 -0.0336537732
## [41] 0.0228291229 0.0524292230 -0.0041648258 -0.0050230176 0.0135754565
## [46] 0.0249270560 -0.0081785632 -0.0183640582 0.1051999776 -0.0229445378
## [51] -0.0228678540 0.0856948453 -0.0607329534 -0.0098096240 0.0548939896
## [56] -0.0254791805 0.0203174878 0.0370291926 0.0397719361 -0.0382685950
## [61] -0.0124796023 0.1097979396 -0.0347700662 -0.0304098636 -0.0477748486
## [66] 0.0193811826 0.0344944596 0.1699930012 -0.0773950092 0.0140066038
## [71] -0.0440997041 -0.0075052207 -0.0346184070 -0.0355643322 0.0001369087
## [76] 0.0870761013 0.0114203479 0.0614805387 -0.0694595747 -0.0823767888
## [81] -0.0749754267 -0.0201832002 0.0065058688 0.0165781899 0.0198502904
## [86] 0.0165448715 0.1479444099 -0.1301835994 -0.0374412980 -0.0489044156
##
## $wnv
## [1] 0.00249719
##
## $se
## [1] 0.04997189 0.05962512 0.05976305 0.06272972 0.06357932 0.06364458
## [7] 0.06404001 0.06411843 0.06413654 0.06418891
##
## $psi
## [1] -0.65089333 -0.08120699 0.38147303 -0.20729935 -0.05766419 0.14219228
## [7] -0.06343914 -0.03049638 0.05187835 -0.01837059forNS = fore.arma.wge(p.ns,phi=estp.ns$phi, n.ahead=10,lastn=TRUE,plot=TRUE,limits=FALSE)
Here we can see the forecast is going to ascillate back towards the mean because of the AR2 model, doing the lookback window, we can see that that this doesn’t seem to align with the realization we have.
Let’s confirm for sure by checking the ASE.
ASEaNS = mean((p.ns[(90-20+1):90] - forNS$f)^2)
ASEaNS # 0.0054## [1] 0.005358956trainingSize = 70
horizon = 10
ASEHolder = numeric()
for( i in 1:(81-(trainingSize + horizon) + 1))
{
forecasts = fore.aruma.wge(p.ns[i:(i+(trainingSize-1))],phi = estp.ns$phi, theta = estp.ns$theta, s = 7, d = 0,n.ahead = horizon)
ASE = mean((p.ns[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
ASEHolder[i] = ASE
}
ASEHolder## [1] 0.01165421 0.01169001hist(ASEHolder)
WindowedASE = mean(ASEHolder)
summary(ASEHolder)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.01165 0.01166 0.01167 0.01167 0.01168 0.01169WindowedASE## [1] 0.01167211Let’s now look to account for the possible seasonality as well as a trend to see if we can get a better ASE.
# Beginning the process of building a seasonal model (s=7)
y=artrans.wge(posperas, phi.tr=c(0,0,0,0,0,0,1))
y_difTwice=artrans.wge(p.ns, phi.tr=c(0,0,0,0,0,0,1))
We can check the white noise with the ACF on trend and seasonal differenced data.
acf(y_difTwice)
ljung.wge(y_difTwice)## Obs -0.5086921 -0.1082985 0.2776785 -0.2555193 0.0553681 0.2680201 -0.4014628 0.2043357 -0.06299046 0.1172208 0.04140118 -0.2462231 0.1937597 -0.09924006 -0.07003595 0.2608693 -0.2185733 -0.0554479 0.2442204 -0.1992782 0.1296738 -0.03764934 -0.1045097 0.09942883## $test
## [1] "Ljung-Box test"
##
## $K
## [1] 24
##
## $chi.square
## [1] 103.0335
##
## $df
## [1] 24
##
## $pval
## [1] 9.087064e-12ljung.wge(y_difTwice, K = 48)## Obs -0.5086921 -0.1082985 0.2776785 -0.2555193 0.0553681 0.2680201 -0.4014628 0.2043357 -0.06299046 0.1172208 0.04140118 -0.2462231 0.1937597 -0.09924006 -0.07003595 0.2608693 -0.2185733 -0.0554479 0.2442204 -0.1992782 0.1296738 -0.03764934 -0.1045097 0.09942883 -0.02804324 0.009471847 0.07446027 -0.1566096 0.07507464 0.0745692 -0.11788 0.06207474 0.02179908 -0.1185527 0.08455902 0.05604227 -0.1051176 0.04851398 -0.0001723275 -0.04692321 0.06820976 -0.01528358 -0.04306204 0.03018357 -0.02891664 0.08508824 -0.03968639 -0.04787157## $test
## [1] "Ljung-Box test"
##
## $K
## [1] 48
##
## $chi.square
## [1] 120.5375
##
## $df
## [1] 48
##
## $pval
## [1] 3.608978e-08Here we can see that we have residuals exceeding our confidence limits and our ljung box test rejects the hypothesis that our model is just white noise.
Now let’s check the AIC of the trend and seasonal models.
aic5.wge(y_difTwice, p = 0:12, type='bic') #selects a aruma(11,1,0)s=7## ---------WORKING... PLEASE WAIT...
##
##
## Error in aic calculation at 10 2
## Five Smallest Values of bic## p q bic
## 34 11 0 -5.549322
## 36 11 2 -5.521321
## 28 9 0 -5.504017
## 35 11 1 -5.498176
## 37 12 0 -5.495912The bic, which punishes model complexity gave us a preferred ARUMA (11,1,0) s= 7. Now let’s introduce the estimate model to get the phis to forecast against.
esty = est.arma.wge(y_difTwice,p=11, q=0)##
## Coefficients of Original polynomial:
## -0.9316 -0.7763 -0.4537 -0.5642 -0.5467 -0.2210 -0.4410 -0.2350 0.0190 0.4947 0.3624
##
## Factor Roots Abs Recip System Freq
## 1-0.5077B+0.9517B^2 0.2668+-0.9898i 0.9755 0.2081
## 1-1.4896B+0.9438B^2 0.7892+-0.6609i 0.9715 0.1110
## 1+1.3305B+0.9173B^2 -0.7252+-0.7511i 0.9578 0.3722
## 1+0.7449B+0.8840B^2 -0.4213+-0.9766i 0.9402 0.3148
## 1+0.8814B -1.1346 0.8814 0.5000
## 1-0.7654B 1.3065 0.7654 0.0000
## 1+0.7376B -1.3557 0.7376 0.5000
##
## esty$phi## [1] -0.93160508 -0.77629346 -0.45374076 -0.56419963 -0.54672712 -0.22102422
## [7] -0.44101332 -0.23499070 0.01900984 0.49472099 0.36241750esty$theta## [1] 0factor.wge(phi=esty$phi)##
## Coefficients of Original polynomial:
## -0.9316 -0.7763 -0.4537 -0.5642 -0.5467 -0.2210 -0.4410 -0.2350 0.0190 0.4947 0.3624
##
## Factor Roots Abs Recip System Freq
## 1-0.5077B+0.9517B^2 0.2668+-0.9898i 0.9755 0.2081
## 1-1.4896B+0.9438B^2 0.7892+-0.6609i 0.9715 0.1110
## 1+1.3305B+0.9173B^2 -0.7252+-0.7511i 0.9578 0.3722
## 1+0.7449B+0.8840B^2 -0.4213+-0.9766i 0.9402 0.3148
## 1+0.8814B -1.1346 0.8814 0.5000
## 1-0.7654B 1.3065 0.7654 0.0000
## 1+0.7376B -1.3557 0.7376 0.5000
##
## The preferred estimated model is (1 + 0.93B + 0.78B^2 + 0.45B^3 + 0.56B^4 + 0.55B^5 + 0.22B^6 + 0.44B^7 + 0.23B^8 - 0.02B^9 - 0.49B^10 - 0.36B^11)(1 - B^7)(1 - B) = a_t
fore.aruma.wge(y_difTwice,phi=esty$phi, theta=esty$theta, d=1,s=7,n.ahead=20, lastn=FALSE, plot=TRUE)
## $f
## [1] -0.117982177 0.078933948 -0.024670862 0.087408015 0.020904860
## [6] 0.029911221 0.046371654 -0.100961496 -0.056119951 0.079508564
## [11] 0.078390927 -0.076578970 0.094271282 0.004249576 0.001067795
## [16] 0.052371081 -0.054422953 0.089037175 -0.075872605 0.053578832
##
## $ll
## [1] -0.3291555 -0.1327327 -0.2390170 -0.1389213 -0.2055821 -0.1978124
## [7] -0.1981435 -0.3982141 -0.3657227 -0.2566606 -0.2938486 -0.4527361
## [13] -0.2842572 -0.3907782 -0.4496703 -0.4115868 -0.5597535 -0.4475741
## [19] -0.6239520 -0.5175877
##
## $ul
## [1] 0.09319118 0.29060064 0.18967524 0.31373733 0.24739186 0.25763488
## [7] 0.29088680 0.19629115 0.25348276 0.41567773 0.45063043 0.29957819
## [13] 0.47279979 0.39927735 0.45180592 0.51632899 0.45090758 0.62564845
## [19] 0.47220675 0.62474541
##
## $resid
## [1] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
## [6] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
## [11] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
## [16] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0681430774
## [21] -0.0706129159 -0.0225526374 -0.0274188020 0.0356701158 0.0332258995
## [26] -0.0011299545 -0.0086625245 -0.0559412915 0.0458540778 0.0668773927
## [31] -0.0331571912 -0.0216007497 -0.0062198633 0.0569491009 0.0523944456
## [36] -0.0935384744 -0.0320721403 0.0413442280 0.0294130617 -0.0081301100
## [41] -0.0832824272 0.0947978243 -0.0403392237 -0.0526683321 0.0939437766
## [46] -0.1079804120 0.0569360292 0.0789274147 -0.2586533369 0.1852108042
## [51] 0.1128448642 -0.1715995049 0.0598045492 0.0003855599 0.0631412755
## [56] 0.0257384618 -0.1125492461 -0.1053374184 0.1338773664 0.0975051230
## [61] -0.0009518359 -0.3156554663 0.2634177102 -0.0179090914 0.0860108224
## [66] -0.1611953978 -0.0633765116 -0.0488457602 0.2972979582 -0.2188026752
## [71] 0.1037161444 -0.0353242712 0.0382993894 -0.0504722564 0.0193915529
## [76] -0.0136427102 0.0686780672 -0.0751314076 -0.0356386053 0.2452012241
## [81] -0.1001504273 -0.0111726594 -0.0538298198
##
## $wnv
## [1] 0.01160823
##
## $se
## [1] 0.1077415 0.1079932 0.1093603 0.1154741 0.1155546 0.1161855 0.1247526
## [8] 0.1516595 0.1579606 0.1715149 0.1899181 0.1919169 0.1931268 0.2015448
## [15] 0.2299684 0.2367132 0.2578217 0.2737813 0.2796323 0.2914115
##
## $psi
## [1] 0.06839492 0.15998949 0.34411769 -0.04001378 0.11223147 0.42168321
## [7] 0.80043719 0.40995670 0.62025073 0.75696874 0.25641297 0.20032771
## [13] 0.53498204 1.02790720 0.52072870 0.94829606 0.85491138 0.52815237
## [19] 0.76121876 0.76716125
##
## $ptot
## [1] 19
##
## $phitot
## [1] 0.06839492 0.15531162 0.32255270 -0.11045887 0.01747252 0.32570289
## [7] 0.78001091 0.13762770 0.09868892 0.15315846 -0.02184462 -0.37989002
## [13] -0.32570289 0.21998909 -0.20602262 -0.25400054 -0.47571115 0.13230349
## [19] 0.36241750forS7 = fore.aruma.wge(y_difTwice,phi=esty$phi, theta=esty$theta, d=1,s=7,n.ahead=20, lastn=TRUE, plot=TRUE)
The seasonal and trend forecast does seem to pass the visual test with our cycling being captured and the trend looking like it’s starting to dip back down.
Let’s check for certain with the ASE.
ASEas7 = mean((y_difTwice[(83-20+1):83] - forS7$f)^2)
ASEas7## [1] 0.01183746trainingSize = 70
horizon = 10
ASEHolder = numeric()
for( i in 1:(81-(trainingSize + horizon) + 1))
{
forecasts = fore.aruma.wge(y_difTwice[i:(i+(trainingSize-1))],phi = esty$phi, theta = esty$theta, s = 7, d = 1,n.ahead = horizon)
ASE = mean((y_difTwice[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
ASEHolder[i] = ASE
}
ASEHolder## [1] 0.01233464 0.01045517hist(ASEHolder)
WindowedASE = mean(ASEHolder)
summary(ASEHolder)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.01046 0.01093 0.01139 0.01139 0.01186 0.01233WindowedASE## [1] 0.0113949Now lets look at all the ASEs from our stationary, trend and seasonal model to see which one performed the best.
print('So, In Summary, the ASE Values Are:')## [1] "So, In Summary, the ASE Values Are:"print(paste('Arma(1,1):', ASEa22))## [1] "Arma(1,1): 0.00258659893879501"print(paste('Arima(2,1,0):', ASEaNS))## [1] "Arima(2,1,0): 0.0053589563896785"print(paste('Aruma(11,1,0)s=7:', ASEas7))## [1] "Aruma(11,1,0)s=7: 0.0118374559290337"Here it appears that our staionary model had the lowest ASE for the AZ data with 0.0026.
Next, we will look into including multiple variables to see hwo that affects positive percentage and the ASE score we achieved.
AZ Var Forecast
For our multivariate forecasts, we will use the VAR model to predict positive percentage in the US and AZ data.
To get started, let’s look into seeing if there’s a lag between other variables and the ccf.
posperas = (as2$positiveIncrease[1:140]/as2$totalTestResultsIncrease[1:140])
posperas[is.na(posperas)] <- 0
as2[is.na(as2)] <- 0
ccf(posperas, as2$death) # lag is -18
ccf(posperas, as2$hospitalizedCurrently) # lag is 0
ccf(posperas, as2$inIcuCurrently) #lag is -15
ccf(posperas, as2$onVentilatorCurrently) # lag is -15
ccf(posperas, as2$deathIncrease) #lag is 12
ccf(posperas, as2$totalTestResultsIncrease) #lag is -1
ccf(posperas, as2$positiveIncrease)
Var Forecast of Stationary data.
Looking at the CCFs of positive percentage with other key variables, we can see that almost all of them have some sort of lagged correlation with the residual variables. So when we run a multivariate forecast, we can account for these lags to strengthen the forecast.
as2[is.na(as2)] <- 0For a cleaning standpoint and to allow our models to run, we will make sure that our NAs are converted to 0s.
library(vars)
as2$posper = posperas
#us2$posper[1:172], us2$death[1:172], us2$hospitalizedCurrently[1:172], us2$inIcuCurrently[1:172], us2$onVentilatorCurrently[1:172], us2$deathIncrease[1:172], us2$totalTestResultsIncrease[1:172], us2$negativeIncrease[1:172], us2$positiveIncrease[1:172]
BSVar1 = VARselect(cbind(as2$posper[1:140], as2$totalTestResultsIncrease[1:140], as2$positiveIncrease[1:140]), type = 'both', lag.max = 5)
BSVar1 # AIC was 2.38e01 for p =3## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 3 3 1 3
##
## $criteria
## 1 2 3 4 5
## AIC(n) 2.404498e+01 2.389207e+01 2.379603e+01 2.389824e+01 2.386973e+01
## HQ(n) 2.417616e+01 2.410196e+01 2.408463e+01 2.426555e+01 2.431575e+01
## SC(n) 2.436778e+01 2.440856e+01 2.450621e+01 2.480211e+01 2.496728e+01
## FPE(n) 2.771049e+10 2.378895e+10 2.162499e+10 2.398008e+10 2.334776e+10Running a var that accounts for the lags between positive percentage, deaths, numbers of people in the ICU, the number of hospitalized patients, number of people on ventilators, number of positives, number of negatives and total test results, we can see that accounting for all the lags the preferred AIC is a p of 3 at an AIC of 3.79e+01.
Now let’s look to see if we can predict the forecast on the var.
#ustest <- data.frame(posper, us2$death, us2$hospitalizedCurrently, us2$inIcuCurrently, us2$onVentilatorCurrently, us2$totalTestResultsIncrease, us2$negativeIncrease, us2$positiveIncrease)
#ustest = cbind(us2$posper[40:140], us2$death[40:140], us2$hospitalizedCurrently[40:140], us2$inIcuCurrently[40:140], us2$onVentilatorCurrently[40:140], us2$totalTestResultsIncrease[40:140], us2$negativeIncrease[40:140], us2$positiveIncrease[40:140])
PosperVAR = VAR(cbind(as2$posper[1:140], as2$totalTestResultsIncrease[1:140], as2$positiveIncrease[1:140]), type = "const", p = 3)## Warning in VAR(cbind(as2$posper[1:140], as2$totalTestResultsIncrease[1:140], : No column names supplied in y, using: y1, y2, y3 , instead.predsar10=predict(PosperVAR,n.ahead=10)
predsar20=predict(PosperVAR,n.ahead=20)To run our VAR model, we will need to ensure that our test results and positive increase variables are both working off of a 140 row dataframe using the p of 3 to run our var against. Next we will predict against a 10 day and 20 day horizon as our short term and long term forecast.
startingPoints10 = us2$posper[131:140]
posperForecasts10 = predsar10$fcst$y1[,1:4] + startingPoints10
startingPoints20 = us2$posper[121:140]
posperForecasts20 = predsar20$fcst$y1[,1:4] + startingPoints20Here we will create two separate predictor variables against our VAR stationary model for the 10 and 20 day horizon.
#dev.off()
library(RColorBrewer)
#fanchart(preds, colors = brewer.pal(n = 8, name= 'Blues'))
plot(seq(1,140,1), as2$posper, type = "l",xlim = c(0,145), ylim = c(0,0.6), ylab = "Positive Percentage", main = "10 Day Positive Percentage Forecast")
lines(seq(131,140,1), as.data.frame(posperForecasts10)$fcst, type = "l", col = "red")
plot(seq(1,140,1), as2$posper, type = "l",xlim = c(0,145), ylim = c(0,0.6), ylab = "Positive Percentage", main = "20 Day Positive Percentage Forecast")
lines(seq(121,140,1), as.data.frame(posperForecasts20)$fcst, type = "l", col = "red")
Using the Var forecast including death, hospitalized patients, patients on ventilators and patients in ICU, we can see that the prediction over Seemed to expect the spikes from 100 to 125 to continue, which it did not in actuality. THe forecasts for the short term and long term trends both seemed to be problematic in the same way.
ASEavarstat = mean((as2$posper[131:140] - posperForecasts10[,1])^2)
ASEavarstat #0.018## [1] 0.01823804ASEavarstat20 = mean((as2$posper[121:140] - posperForecasts20[,1])^2)
ASEavarstat20 #0.022## [1] 0.02191405Now lets do a var assuming a trend model (which ended up being the best in a univariate study)
Var forecast with trend
pospertrend = artrans.wge(as2$posper, 1)
testresultstrends = artrans.wge(as2$totalTestResultsIncrease, 1)
positivetesttrend = artrans.wge(as2$positiveIncrease, 1)
Here, we will need to difference the two variables we included in the stationary VAR model as well as the positive percentage.
Now that we have differenced variables, lets see which gives us the best AIC score
BSVar1 = VARselect(cbind(pospertrend[1:139], testresultstrends[1:139], positivetesttrend[1:139]), type = 'both', lag.max = 10)
BSVar1 # AIC was 2.39e01 for p =9## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 9 2 2 9
##
## $criteria
## 1 2 3 4 5
## AIC(n) 2.445036e+01 2.414332e+01 2.424580e+01 2.412701e+01 2.417328e+01
## HQ(n) 2.458548e+01 2.435951e+01 2.454305e+01 2.450533e+01 2.463267e+01
## SC(n) 2.478290e+01 2.467538e+01 2.497738e+01 2.505811e+01 2.530390e+01
## FPE(n) 4.156310e+10 3.058583e+10 3.391226e+10 3.015391e+10 3.164682e+10
## 6 7 8 9 10
## AIC(n) 2.399090e+01 2.397968e+01 2.403631e+01 2.385806e+01 2.389174e+01
## HQ(n) 2.453137e+01 2.460121e+01 2.473891e+01 2.464174e+01 2.475649e+01
## SC(n) 2.532105e+01 2.550935e+01 2.576550e+01 2.578677e+01 2.601998e+01
## FPE(n) 2.644856e+10 2.625799e+10 2.793359e+10 2.352901e+10 2.453709e+10Like the earlier model, we end up with an AIC preference of 9.
PosperVAR = VAR(cbind(pospertrend[1:139], testresultstrends[1:139], positivetesttrend[1:139]), type = "const", p = 9)## Warning in VAR(cbind(pospertrend[1:139], testresultstrends[1:139], positivetesttrend[1:139]), : No column names supplied in y, using: y1, y2, y3 , instead.predstrend10=predict(PosperVAR,n.ahead=10)
predstrend20=predict(PosperVAR,n.ahead=20)LIke previously, we will run a VAR in the test results, positive test trend and pospertrend within the observations 1-139. Additionally, we will be looking at a horizon of 10 and 20.
startingPoints10 = pospertrend[130:139]
posperForecasts10 = predstrend10$fcst$y1[,1:3] + startingPoints10
startingPoints20 = pospertrend[120:139]
posperForecasts20 = predstrend20$fcst$y1[,1:3] + startingPoints20This will set up our trend forecasts from the 10 day to 20 day horizon so that we can plot them.
#dev.off()
library(RColorBrewer)
#fanchart(preds, colors = brewer.pal(n = 8, name= 'Blues'))
plot(seq(1,139,1), pospertrend, type = "l",xlim = c(0,145), ylim = c(-0.5,0.5), ylab = "Positive Percentage", main = "10 Day Positive Percentage Forecast")
lines(seq(130,139,1), as.data.frame(posperForecasts10)$fcst, type = "l", col = "red")
plot(seq(1,139,1), pospertrend, type = "l",xlim = c(0,145), ylim = c(-0.5,0.5), ylab = "Positive Percentage", main = "20 Day Positive Percentage Forecast")
lines(seq(120,139,1), as.data.frame(posperForecasts20)$fcst, type = "l", col = "red")
Our VAR forecast using a trend model did a really good job of capturing the seasonality trend in the 10 and 20 day forecast. It seemed to capture perfectly some of the spikes from 135 to 139 and 125 - 130. This was apparent in the 10 and 20 day forecast.
ASEavartrend = mean((pospertrend[130:139] - posperForecasts10[,1])^2)
ASEavartrend #0.0017## [1] 0.001688938ASEavartrend20 = mean((pospertrend[120:139] - posperForecasts20[,1])^2)
ASEavartrend20 #0.0017## [1] 0.001239529We can see with a trend model, we get a better ASE than than the stationary model we ran previously. Both our 10 and 20 day ASEs had strong scores which were better than anything we achieved previously.
Now lets try a seasonal model accounting for 7 days on top of the trend data we previously ran.
pospertrend2 = artrans.wge(pospertrend, c(rep(0,6),1))
testresultstrends2 = artrans.wge(testresultstrends, c(rep(0,6),1))
positivetesttrend2 = artrans.wge(positivetesttrend, c(rep(0,6),1))
Here we add a 7 day seasonality to our trend forecast so that our models take on an Airline like model. Looking at our transformed realization, there doesn’t seem to be a major difference, which suggests that we would probably not see an improved performance to our trend model.
BSVar1 = VARselect(cbind(pospertrend2[1:132], testresultstrends2[1:132], positivetesttrend2[1:132]), type = 'both', lag.max = 10)
BSVar1 # AIC was 2.41e01 for p =9## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 9 9 4 9
##
## $criteria
## 1 2 3 4 5
## AIC(n) 2.598000e+01 2.549624e+01 2.540144e+01 2.496787e+01 2.491157e+01
## HQ(n) 2.612003e+01 2.572028e+01 2.570950e+01 2.535995e+01 2.538767e+01
## SC(n) 2.632476e+01 2.604785e+01 2.615990e+01 2.593319e+01 2.608374e+01
## FPE(n) 1.918806e+11 1.183369e+11 1.077318e+11 6.994088e+10 6.627270e+10
## 6 7 8 9 10
## AIC(n) 2.471481e+01 2.438685e+01 2.429823e+01 2.409346e+01 2.416891e+01
## HQ(n) 2.527493e+01 2.503099e+01 2.502638e+01 2.490563e+01 2.506510e+01
## SC(n) 2.609384e+01 2.597273e+01 2.609096e+01 2.609305e+01 2.637535e+01
## FPE(n) 5.462556e+10 3.953866e+10 3.641083e+10 2.990431e+10 3.256714e+10Like previously, we found that that the AIC gave us a p of 9 with a top AIC score of 2.41.
PosperVAR = VAR(cbind(pospertrend2[1:132], testresultstrends2[1:132], positivetesttrend2[1:132]), type = "const", p = 9)## Warning in VAR(cbind(pospertrend2[1:132], testresultstrends2[1:132], positivetesttrend2[1:132]), : No column names supplied in y, using: y1, y2, y3 , instead.preds10=predict(PosperVAR,n.ahead=10)
preds20=predict(PosperVAR,n.ahead=20)Next we will run our var forecasts against our 3 variables that were trend and seasonally differenced at a p = 9.
startingPoints10 = pospertrend2[123:132]
posperForecasts10 = preds10$fcst$y1[,1:3] + startingPoints10
startingPoints20 = pospertrend2[113:132]
posperForecasts20 = preds20$fcst$y1[,1:3] + startingPoints20Here we will prepare our 10 and 20 day forecast to plot and measure ASE.
#dev.off()
library(RColorBrewer)
#fanchart(preds, colors = brewer.pal(n = 8, name= 'Blues'))
plot(seq(1,132,1), pospertrend2, type = "l",xlim = c(0,140), ylim = c(-0.5,0.5), ylab = "Positive Percentage", main = "10 Day Positive Percentage Forecast")
lines(seq(123,132,1), as.data.frame(posperForecasts10)$fcst, type = "l", col = "red")
plot(seq(1,132,1), pospertrend2, type = "l",xlim = c(0,140), ylim = c(-0.5,0.5), ylab = "Positive Percentage", main = "20 Day Positive Percentage Forecast")
lines(seq(113,132,1), as.data.frame(posperForecasts20)$fcst, type = "l", col = "red")
From our 10 and 20 day forecast, we can see that our seasonality over emphasized the spikes that took place at 110 to build it’s forecast for the days to come after.
ASEavarairplane = mean((pospertrend2[123:132] - posperForecasts10[,1])^2)
ASEavarairplane #0.014## [1] 0.01478145ASEavarairplane20 = mean((pospertrend2[113:132] - posperForecasts20[,1])^2)
ASEavarairplane20 #0.009## [1] 0.009607299In doing so, we can see that it led to one of the poorest ASE scores that we came across.
print('So, In Summary, the ASE Values Are:')## [1] "So, In Summary, the ASE Values Are:"print(paste('Arma(2,2):', ASEa22))## [1] "Arma(2,2): 0.00258659893879501"print(paste('Arima(2,1,0):', ASEaNS))## [1] "Arima(2,1,0): 0.0053589563896785"print(paste('Aruma(9,1,0)s=7:', ASEas7))## [1] "Aruma(9,1,0)s=7: 0.0118374559290337"print(paste('Var Stationary:', ASEavarstat))## [1] "Var Stationary: 0.0182380429222605"print(paste('Var Trend:', ASEavartrend))## [1] "Var Trend: 0.00168893770664753"print(paste('Var Airplane:', ASEavarairplane))## [1] "Var Airplane: 0.0147814492897754"So far we can see that our bsest performing model thus far was a Var model that accounted for a trend in the analysis. This is a bit of a disconnect from our Fore Arma model forecasts that chose the stationary model as the top performer. That said, looking at the ARMA models and the VAR models together, the Trend forecast models had the lowest average.
MLP
Since our trend model was the best overall performer across our ARMA forecasts and Var forecasts, we will do an MLP of our trend forecast.
We will break up our US2 dataset into two separate horizons at 10 days and 20 days and then we will difference such that you can build a trend.
library(nnfor)
library(forecast)
library(vars)
library(tswge)
as2Small = as2[1:130,]
as2Large = as2[1:120,]
pospersmall = artrans.wge(as2Small$posper, phi.tr=(1))
testresultssmall = artrans.wge(as2Small$totalTestResultsIncrease, phi.tr=(1))
positivesmall = artrans.wge(as2Small$positiveIncrease, phi.tr=(1))
posperlarge = artrans.wge(as2Large$posper, phi.tr=(1))
testresultslarge = artrans.wge(as2Large$totalTestResultsIncrease, phi.tr=(1))
positivelarge = artrans.wge(as2Large$positiveIncrease, phi.tr=(1))
as2SmallDF = data.frame(test = ts(testresultssmall, frequency = 1), positive = ts(positivesmall, frequency = 1))
as2LargeDF = data.frame(test = ts(testresultslarge, frequency = 1), positive = ts(positivelarge, frequency = 1))In doing so, we will also difference positive percentage, test results increase and positive increase so that we have a trend. We will separate our variabiles into small and large forecasts.
From here, we will fit it into a neural net model.
fit.mlp.small = mlp(ts(pospersmall, frequency = 1), allow.det.season = FALSE, reps = 20, comb = 'mean', xreg = as2SmallDF)
#fit.mlp.small = mlp(ts(us2Small$deathIncrease, frequency = 1), allow.det.season = FALSE, reps = 20, comb = 'mean', xreg = us2SmallDF)
fit.mlp.small## MLP fit with 5 hidden nodes and 20 repetitions.
## Univariate lags: (1,2)
## Forecast combined using the mean operator.
## MSE: 0.0198.plot(fit.mlp.small)
fit.mlp.large = mlp(ts(posperlarge, frequency = 1), allow.det.season = FALSE, reps = 20, comb = 'mean', xreg = as2LargeDF)
#fit.mlp.large = mlp(ts(us2Large$deathIncrease, frequency = 1), allow.det.season = FALSE, reps = 20, comb = 'mean', xreg = us2LargeDF)
fit.mlp.large## MLP fit with 5 hidden nodes and 20 repetitions.
## Univariate lags: (1,2)
## Forecast combined using the mean operator.
## MSE: 0.0218.plot(fit.mlp.large)
The MLP output gave us 5 notes and 20 repititions, with two regressors included. It did the same for both the short term and long term forecast.
From here we will create a dataframe using the key features of the regressors, and then forecast against it using a horizon of 10 days.
testdif = artrans.wge(as2$totalTestResults, phi.tr=(1))
posdif = artrans.wge(as2$positiveIncrease, phi.tr=(1))
ASDF = data.frame(test = ts(testdif), positive = ts(posdif))
fore.mlp10 = forecast(fit.mlp.small, h = 10, xreg = ASDF)
fore.mlp20 = forecast(fit.mlp.large, h = 20, xreg = ASDF)
plot(fore.mlp10)
plot(fore.mlp20)
Within our 20 repetitions, we can see that our 10 day and 20 day MLP seemed to pass the eye test of continuing the data that appears before it.
Now let’s check the ASE of our short term and long term forecast.
posperdif = artrans.wge(as2$posper, phi.tr=(1))
ASEmlp10 = mean((posperdif[130:139] - fore.mlp10$mean)^2)
ASEmlp10 #0.0077## [1] 0.007599468ASEmlp20 = mean((posperdif[120:139] - fore.mlp20$mean)^2)
ASEmlp20 #0.0076## [1] 0.004901224Our ASE of the MLP model seemed to be comparable to the trend models that we saw previously. As both were in the 0.001 to 0.007 range.
Ensemble Model
Since our trend differenced ARMA and Var models did the best and informed our MLP model, we will use the two to inform our ensemble model to see if it does better than 0.001 in ASE.
ensemble10 = (predstrend10$fcst$y1[,1] + fore.mlp10$mean)/2
ensemble20 = (predstrend20$fcst$y1[,1] + fore.mlp20$mean)/2plot(seq(1,139,1), pospertrend, type = "l",xlim = c(0,145), ylim = c(-0.5,0.5), ylab = "Positive Percentage", main = "10 Day Positive Percentage Forecast Ensemble")
lines(seq(130,139,1), ensemble10, type = "l", col = "red")
plot(seq(1,139,1), pospertrend, type = "l",xlim = c(0,145), ylim = c(-0.5,0.5), ylab = "Positive Percentage", main = "20 Day Positive Percentage Forecast Ensemble")
lines(seq(120,139,1), ensemble20, type = "l", col = "red")
Our ensemble model of the mlp and the var appeared to split the difference a bit between our model that captured the spikes really well to one that showed consistent oscillating trend to zero.
Ensemble ASE
ASE10 = mean((pospertrend[130:139] - ensemble10)^2)
ASE10 # 0.008## [1] 0.008332369ASE20 = mean((pospertrend[120:139] - ensemble20)^2)
ASE20 #0.006## [1] 0.005055811Our Ensemble Model did not improve upon our Var trend model forecast which garnered an ASE of 0.001. Here, we did significantly higher at 0.008 and 0.006 respectively.
print('So, In Summary, the ASE Values Are:')## [1] "So, In Summary, the ASE Values Are:"print(paste('Arma(2,2):', ASEa22))## [1] "Arma(2,2): 0.00258659893879501"print(paste('Arima(2,1,0):', ASEaNS))## [1] "Arima(2,1,0): 0.0053589563896785"print(paste('Aruma(9,1,0)s=7:', ASEas7))## [1] "Aruma(9,1,0)s=7: 0.0118374559290337"print(paste('Var Stationary:', ASEavarstat))## [1] "Var Stationary: 0.0182380429222605"print(paste('Var Trend:', ASEavartrend))## [1] "Var Trend: 0.00168893770664753"print(paste('Var Airplane:', ASEavarairplane))## [1] "Var Airplane: 0.0147814492897754"print(paste('MLP:', ASEmlp10))## [1] "MLP: 0.00759946773558966"print(paste('Ensemble:', ASE10))## [1] "Ensemble: 0.00833236875504159"Overall, our Best performing forecast of the Arizona positive percentage was the Var Trend Model in both the short term and long term, which is a bit more inaccurate than our US forecast. The reason for this could be due to the shape of the positive percentate realization for the US compared to Arizona, where the realizations in the US data has decreased close to zero by the time of the study while the AZ data was still close to it’s apex.